AFRL-ML-WP-TR-2002-4060 


AEROSPACE  COMPOSITE  MATERIALS 
Delivery  Order  0002:  Development  and 
Validation  of  Micromechanical  Models  for 
Composites 


V.  Buryachenko 


E.  larve 


R.Kiin 


S.  Sihn 


G.  P.  Tandon 


University  of  Dayton  Research  Institute 
300  College  Park  Avenue 
Dayton,  OH  45469-0168 


MARCH  2002 


Interim  Report  for  30  November  2000  -  29  November  2001 


Approved  for  public  release;  distribution  is  unlimited. 


MATERIALS  AND  MANUFACTURING  DIRECTORATE 
AIR  FORCE  RESEARCH  LABORATORY 
AIR  FORCE  MATERIEL  COMMAND 
WRIGHT-PATTERSON  AIR  FORCE  BASE,  OH  45433-7750 


NOTICE 


Using  government  drawings,  specifications,  or  other  data  included  in  this  document  for  any 
purpose  other  than  government  procurement  does  not  in  any  way  obligate  the  U.S.  Government. 
The  fact  that  the  government  formulated  or  supplied  the  drawings,  specifications,  or  other  data 
does  not  license  the  holder  or  any  other  person  or  corporation;  or  convey  and  rights  or  permission 
to  manufacture,  use,  or  sell  any  patented  invention  that  may  relate  to  them. 

This  report  has  been  reviewed  by  the  Office  of  Public  Affairs  (ASC/PA)  and  is  releasable  to  the 
National  Technical  Information  Service  (NTIS).  At  NTIS,  it  will  be  available  to  the  general 
public,  including  foreign  nations. 

This  technical  report  has  been  reviewed  and  is  approved  for  publication. 


KARLA  L.  STRONG,  Project  ^gineer 
Structural  Materials  Branch  ^ 
Nonmetallic  Materials  Division 


TIA  BENSON  TOLLE,  Acting  Chief 
Structural  Materials  Branch 
Nonmetallic  Materials  Division 


ROBERT  M.  SUSNIK,  Deputy  Chief 
Nonmetallic  Materials  Division 
Materials  and  Manufacturing  Directorate 


Copies  of  this  report  should  not  be  returned  unless  return  is  required  by  security  considerations, 
contractual  obligations,  or  notice  on  a  specific  document. 


REPORT  DOCUMENTATION  PAGE 


Form  Approved 
0MB  No.  0704-0188 


The  public  reporting  burden  for  this  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources,  searching  existing  data  sources, 
gathering  and  maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information,  including 
suggestions  for  reducing  this  burden,  to  Department  of  Defense,  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports  (07040188),  1215  Jefferson  Davis  Highway,  Suite  1204, 
Arlington,  VA  22202-4302.  Respondents  should  be  aware  that  notwithstanding  any  other  provision  of  law,  no  person  shall  be  subject  to  any  penalty  for  failing  to  comply  with  a  collection  of  information  if  it  does  not 
display  a  currently  valid  0MB  control  number.  PLEASE  DO  NOT  RETURN  YOUR  FORM  TO  THE  ABOVE  ADDRESS. 


1.  REPORT  DATE  (DD-MM-YY)  2.  REPORT  TYPE  3.  DATES  COVERED  (From  -  To) 

March  2002  Interim  1 1/30/2000  -  1 1/29/2001 


5a.  CONTRACT  NUMBER 
F33615-00-D-5006 


5b.  GRANT  NUMBER 


5c.  PROGRAM  ELEMENT  NUMBER 
62102F 


5d.  PROJECT  NUMBER 

4347 


5e.  TASK  NUMBER 

32 _ 

5f.  WORK  UNIT  NUMBER 
02 

8.  PERFORMING  ORGANIZATION 
REPORT  NUMBER 

UDR-TR-2002-00033 


4.  TITLE  AND  SUBTITLE 

AEROSPACE  COMPOSITE  MATERIALS 

Delivery  Order  0002:  Development  and  Validation  of  Micromechanical  Models  for 
Composites 


6.  AUTHOR(S) 

V.  Buryachenko 
E.  larve 

R.  Kim 

S.  Sihn 

G.  P.  Tandon 

7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

University  of  Dayton  Research  Institute 
300  College  Park  Avenue 
Dayton,  Ohio  45469-0168 


9.  SPONSORING/MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES)  10.  SPONSORING/MONITORING 

^  AGENCY  ACRONYM(S) 

Materials  and  Manutacturmg  Directorate 

Air  Force  Research  Laboratory  AFRL/MLBCO _ 

Air  Force  Materiel  Command  11.  SPONSORING/MONITORING 

Wright-Patterson  Air  Force  Base,  OH  45433-7750  AGENCY  REPORT  NUMBER(S) 

AFRL-ML-WP-TR- 2002-4060 

12.  DISTRIBUTION/AVAILABILITY  STATEMENT 
Approved  for  public  release;  distribution  is  unlimited. 


13.  SUPPLEMENTARY  NOTES 
Report  contains  color. 

14.  ABSTRACT  (Maximum  200  Words) 

Mechanics  modeling  and  verification  of  a  variety  of  systems  continues  to  be  examined.  Systems  successfully  examined  include  the 
single  filament  cruciform  test,  identification  of  failure  modes  in  composites,  woven  fabric  composite,  carbon  foams,  the  fiber 
pushout  test,  and  slanted  free-edge  plies.  Techniques  used  include  concentric  cylinder  modeling,  finite  element  method,  B -spline 
approximation,  and  in  situ  observation  of  damage  during  testing  via  scanning  electron  microscopy  (SEM).  The  iteration  method  was 
found  to  be  superior  to  the  Eourier  transform  method  of  modeling  inhomogeneous  composites  except  where  an  analytical  solution  is 
required.  This  knowledge  was  then  applied  to  functionally  graded  materials  and  to  computer-generated  models  of  locally 
nonuniform  composites.  The  radial  distribution  function  for  fiber  locations  and  second-order  intensity  functions  are  needed  to 
predict  properties  beyond  estimating  their  bounds. 

15.  SUBJECT  TERMS 

mechanics  modeling,  micromechanics,  elastic  modulus,  tensile  strength.  Young’s  modulus,  cruciform  test,  epoxy  composites, 
carbon  foams,  axisymmetric  damage  model,  finite  element  method,  EEM,  B-spline  approximation,  B-SAM,  woven  composite, 
fiber  pushout  test,  free-edge  analysis,  micro  structure,  SEM,  random  structures,  iteration  method,  Eourier  transform  method, 
radial  distribution  function,  statistical  second -order  quantities 


16.  SECURITY  CLASSIFICATION  OF: 

a.  REPORT  b.  ABSTRACT  c.  THIS  PAGE 

Unclassified  Unclassified  Unclassified 


17.  LIMITATION 

18.  NUMBER  OF 

OF  ABSTRACT: 

PAGES 

SAR 

170 

19a.  NAME  OF  RESPONSIBLE  PERSON  (Monitor) 
Karla  L.  Strong 

19b.  TELEPHONE  NUMBER  (Include  Area  Code) 
(937)  255-3104 


Standard  Form  298  (Rev.  8-98) 

Prescribed  by  ANSI  Std.  Z39-I8 


1 


TABLE  OF  CONTENTS 


Section  Page 


EXECUTIVE  SUMMARY  1 

1  MIXED-MODE  FAILURE  CRITERIA  USING  CRUCIFORM 

GEOMETRY  4 

I .  I  SUMMARY  AND  RECOMMENDATIONS  1 2 

2  DAMAGE  EVOLUTION  AND  FAILURE  MODELING  IN 

UNIDIRECTIONAL  GRAPHITE/BMI  COMPOSITE  13 

2.1  SUMMARY  AND  RECOMMENDATIONS  23 

3  MICROMECHANICAL  RESPONSE  OF  A  MULTIPHASE 

COMPOSITE  24 

4  SINGLE-FIBER  PUSHOUT  CONSIDERING  INTERFACIAL 

FRICTION  AND  ADHESION  3 1 

5  PARAMETRIC  DESIGN,  ANALYSIS,  AND  TESTING  OF  A 
MODEL  WOVEN  COMPOSITE  WITH  A  CRUCIFORM 

GEOMETRIC  CONFIGURATION  33 

5.1  NUMERICAL  RESULTS  AND  PARAMETRIC  STUDY  35 

5.1.1  Unit  Cell  of  Plain-Weave  Composite  35 

5.1.2  One-Dimensional  Model  Laminate  with  Cruciform 

Geometry  38 

5.2  EXPERIMENT  AND  VALIDATION  42 

5.3  SUMMARY  AND  RECOMMENDATIONS  44 

6  ASYMPTOTIC  STRESS  ANALYSIS  OF  LAMINATED 

COMPOSITES  WITH  SLANTED  FREE  EDGE  47 

6.1  ASYMPTOTIC  ANALYSIS  48 

6.2  FULL-FIELD  NUMERICAL  CALCULATION  55 

6.3  RESULTS  59 

6.4  SUMMARY  AND  RECOMMENDATIONS  77 

7  MECHANICAL  MODELING  AND  PREDICTION  OF  BULK 

PROPERTIES  OF  OPEN-CELL  CARBON  FOAM  79 

7. 1  MODELING  OF  CARBON  FOAM  8 1 


iii 


TABLE  OF  CONTENTS  (Continued) 


Section  Page 


7.1.1  Characterization  of  Geometry  of  Unit  Cell  of 

Carbon  Foam  81 

7.1.2  Determination  o  f  B  oundary  Condition  87 

7.2  NUMERICAL  CALCULATION  88 

7.2. 1  Prediction  Method  of  Effective  Bulk  Moduli  of  Foam  88 

7.2.2  Results:  Polyurethane  Foam  94 

7.2.3  Results:  Carbon  Foam  with  Isotropic  Properties  95 

7.2.4  Results:  Carbon  Foam  with  Anisotropic  Properties  97 

7.3  SUMMARY  AND  RECOMMENDATIONS  102 

8  SCANNING  ELECTRON  MICROSCOPY  INVESTIGATION  OF 
DEFORMATION  IN  THE  VICINITY  OF  A  MODE  II  CRACK  TIP 

IN  UNIDIRECTIONAL  TOUGHENED  EPOXY  COMPOSITES  104 

8.1  SPECIMEN  PREPARATION  104 

8.2  EXPERIMENTAL  PROCEDURE  105 

8.3  RESULTS  AND  DISCUSSION  105 

8.4  CONCLUSIONS  109 

9  NONLOCAL  MODELS  OF  STRESS  FIELD  CONCENTRATIONS 
AND  EFFECTIVE  THERMOELASTIC  PROPERTIES  OF 

RANDOM  STRUCTURE  COMPOSITES  [48]  110 

9. 1  THE  GENERAL  SCHEME  1 1 3 

9.2  NUMERICAL  RESULTS  115 

9.3  CONCLUDING  REMARKS  119 

10  MICROMECHANICS  OF  NONLOCAL  EFFECTS  IN  GRADED 

MATERIALS  122 

10. 1  GENERAL  SCHEME  OF  THE  APPROACH  PROPOSED  125 

10.2  PREVIOUS  NUMERICAL  RESULTS  125 

1 1  QUANTITATIVE  DESCRIPTION  AND  NUMERICAL  SIMULATION 
OF  RANDOM  MICROSTRUCTURES  OF  COMPOSITES  AND 

THEIR  EFFECTIVE  ELASTIC  MODULI  [59 [  1 27 

11.1  STATISTICAL  DESCRIPTION  AND  NUMERICAL 

SIMULATION  OF  RANDOM  STRUCTURE  COMPOSITES  1 30 

11.2  EFFECTIVE  ELASTIC  PROPERTIES  140 

11.3  SUMMARY  AND  RECOMMENDATIONS  141 


IV 


TABLE  OF  CONTENTS  (Concluded) 


Section  Page 

12  PUBLICATIONS  AND  PRESENTATIONS  144 

REFERENCES  147 

LIST  OF  ACRONYMS  152 


V 


LIST  OF  FIGURES 

Figure  Page 

1  Cruciform  Specimen  with  Fiber  Oriented  at  Angle  0  to  the  Loading 

Direction  4 

2  (a)  Visual  Observation  of  Interfacial  Debonding  in  a  45 -degree  Off- 
Axis  Specimen;  (b)  Magnified  View  of  Illuminated  White  Region 

Shown  in  (a)  6 

3  Variation  of  Interfacial  Radial  Stress  Concentration  Factor  with  Off- 

Axis  Angle  8 

4  Variation  of  Interfacial  Shear  Stress  Concentration  Factor  with  Off- 

Axis  Angle  8 

5  Failure  Envelope  for  Off-Axis  Cruciform  Specimens  1 1 

6  Composite  Stress- Strain  Curve  for  IM7/5250-4  Composite  16 

7  Acoustic  Emission  Count  as  a  Eunction  of  Load  Level  16 

8  Typical  Fiber  Breaks  Observed  During  Tensile  Loading  of  IM7/5250-4 

Composite  16 

9  Large  Row  of  Fiber  Breaks  in  a  Single  Plane  Observed  Close  to 

Failure  Load  17 

10  Histogram  of  Measured  Distanee  between  Multiple  Fiber  Breaks  17 

1 1  Matrix  Cracks  Bridged  by  Fibers  17 

12  Interfacial  Debonding  Accompanied  by  Fiber  Break  18 

13  Measured  Distribution  of  IM7  Fiber  Tensile  Strength  20 

14  Shear  Stress  at  the  Interface  22 

15  Fiber  Axial  Stress  at  the  Interfaee  22 

16  Radial  Stress  at  the  Interface  23 

17  Matrix  Axial  Stress  at  the  Interfaee  23 

18  Microstructure  of  Woven  Nextel  720/AS  27 


VI 


LIST  OF  FIGURES  (Continued) 

Figure  Page 

19  Magnified  View  of  Microstructure  Showing  Matrix  Cracks  Due  to 

Processing  27 

20  Effective  Thermoelastic  Moduli  of  Woven  Nextel  720/AS  Composite 

as  a  Function  of  Transverse  Crack  Spacing  27 

21  Overlapping  Cubic  B-Spline  Functions  with  Seven  Subdivisions  34 

22  Unit  Cell  of  a  Plain-Weave  Composite  36 

23  Interfacial  Normal  and  Shear  Stress  Distributions  of  a  Unit  Cell  of  a 

Plain- Weave  Textile  Composite  37 

24  Dimensions  and  Geometry  of  One  Quarter  of  a  Cruciform  Specimen  39 

25  Normalized  In-Plane  Longitudinal  and  Out-of-Plane  Shear  Stress 
Concentration  Factors  Along  the  Middle  of  Transverse  Fill  Yarn  in 

1-D  Model  Laminate  40 

26  Ratio  of  Stress  Concentration  Factors  of  Woven  to  those  of  Fillet 

with  Various  Dimensions  of  the  Cruciform  Specimen  Geometry  41 

27  Ratio  of  Stress  Concentration  Factors  of  Woven  to  those  of  Fillet 

with  Various  Yarn  Waviness  Ratio  43 

28  Comparison  of  Measured  and  Modeled  Cross-Section  Profiles  of 

Model  Woven  Composite  44 

29  Axial  Strain  Comparison  between  B-SAM  Predictions  and  Moire 

Interferometry  Measurements  45 

30  Global  and  Local  Coordinate  Systems  in  Asymptotic  Analysis  49 

31  Laminated  Composite  Plate  with  a  Slanted  Edge  56 

32  Overlapping  Cubic  B-Spline  Functions  with  Seven  Subdivisions  57 

33  Mapping  of  Homogeneous  Domain  in  x-Space  to  a  Unit  Cube  in 

^-Domain  57 

34  Distribution  of  Real  and  Imaginary  Parts  of  the  Roots  with  Various 

Slanted  Angles  60 


vii 


LIST  OF  FIGURES  (Continued) 

Figure 

Page 

35 

Distribution  of  Magnitudes  of  Singular  Stress  Components  with  Two 
Roots  for  (|)  =  1° 

62 

36 

Stress  Distribution  in  Global  xyz-Axes  by  Asymptotic  and  Numerical 
Analyses  near  Singular  Point  Along  Straight  Lines  Parallel  to  the 
Interface  between  Two  Plies 

63 

37 

Stress  Distribution  in  Local  rst- Axes  with  \|/  =  1°  by  Asymptotic  and 
Numerical  Analyses  Along  Contours  of  Constant  Radii  with  the 

Center  at  the  Singular  Point 

64 

38 

Distribution  of  Magnitudes  of  Singular  Stress  Components  with  a 
Complex  Root,  A,i  =  0.9646  +  0.0968  i,  for  (])  =  60° 

67 

39 

Stress  Distribution  in  Gbbal  xyz-Axes  by  Asymptotic  and  Numerical 
Analyses  Near  Singular  Point  Along  Straight  Lines  Parallel  to  the 
Interface  between  Two  Plies 

68 

40 

Stress  Distribution  in  Local  rst- Axes  with  (|)  =  30°  by  Asymptotic  and 
Numerical  Analyses  Along  Contours  of  Constant  Radii  with  the 

Center  at  the  Singular  Point 

69 

41 

Distribution  of  Magnitudes  of  Singular  Stress  Components  with  Two 
Roots  for  (|)  =  45° 

70 

42 

Stress  Distribution  in  Global  xyz-Axes  by  Asymptotic  and  Numerical 
Analyses  Near  Singular  Point  Along  Straight  Lines  Parallel  to  the 
Interface  between  Two  Plies 

71 

43 

Stress  Distribution  in  Local  rst- Axes  with  \|/  =  45°  by  Asymptotic  and 
Numerical  Analyses  Along  Contours  of  Constant  Radii  with  the 

Center  at  the  Singular  Point 

72 

44 

Distribution  of  Magnitudes  of  Singular  Stress  Components  with  Two 
Roots  for  (|)  =  60° 

74 

45 

Stress  Distribution  in  Global  xyz-Axes  by  Asymptotic  and  Numerical 
Analyses  Near  Singular  Point  Along  Straight  Lines  Parallel  to  the 
Interface  between  Two  Plies 

75 

viii 

LIST  OF  FIGURES  (Continued) 


Figure  Page 

46  Stress  Distribution  in  Local  rst- Axes  with  (|)  =  60°  by  Asymptotic  and 
Numerical  Analyses  Along  Contours  of  Constant  Radii  with  the 

Center  at  the  Singular  Point  76 

47  Microstructure  of  Carbon  Foam  80 

48  Points  and  Lines  in  a  Cube  and  a  Tetrahedron  for  Generating  a  Unit 

Cell  of  Carbon  Foam  82 

49  A  Tetrahedron  and  Spheres  to  Generate  a  Unit  Cell  of  Carbon  Foam  83 

50  Unit  Cells  of  Carbon  Foam  with  Various  Porosities  83 


5 1  One  of  Four  Subtetrahedra  Considering  Symmetry  in  Figure  48  84 

52  Cross-Sectional  Slices  During  Bubble-Forming  Process  (R  is  the 

Radius  of  Spheres)  84 

2V6 

53  Intersection  of  a  Subtetrahedron  and  Bubble  Spheres  ( 4^  <R<  a)  85 

54  A  Sub- subtetrahedron  Discretized  into  Finite  Volumes  Along  a 

Ligament  Axis  86 

55  3-D  Mesh  versus  1-D  Beam  Mesh  of  the  Unit  Cell  of  Carbon  Foam  87 


56  An  Imaginary  Tetrahedron  inside  which  a  Unit  Cell  of  the  Foam  is 

Located 


89 


57  Relation  between  Bulk  Density,  Porosity  and  Bubble  Size  for 

Open- Cell  Foams 


90 


58  Coordinate  Systems  for  Rotating  Unit  Cell  for  Simulating  Random 

Orientation 


92 


59  Variation  of  Young’s  Moduli  of  Unit  Cells  of  Isotropic  Foam  on  the 
Changes  of  Angles  in  Cartesian  Coordinate  Axes,  with  the  First 

Rotation  in  the  y-Axis  and  the  Second  Rotation  in  the  z-Axis  92 

60  Variation  of  Young’s  Moduli  of  Unit  Cells  of  Isotropic  Foam  on  the 
Changes  of  Angles  in  Cartesian  Coordinate  Axes,  with  the  First 

Rotation  in  the  z-Axis  and  the  Second  Rotation  in  the  y-Axis  92 


IX 


LIST  OF  FIGURES  (Continued) 

Figure  Page 

61  Variation  of  Young’s  Moduli  of  Unit  Cells  of  Isotropic  Foam  on  the 

Change  of  Angles  in  Spherical  Coordinate  Axes  93 

62  Prediction  of  Bulk  Young’s  Moduli  of  Polyurethane  Foam  at 

Various  Densities  95 

63  Prediction  of  Bulk  Poisson’s  Ratio  of  Polyurethane  Foam  at  Various 

Densities  96 

64  Prediction  of  Bulk  Tensile  Moduli  of  Carbon  Foam  at  Various  Densities  97 

65  Prediction  of  Bulk  Tensile  Moduli  of  Carbon  Foam  at  Various  Densities 

with  Several  Assumed  Ligament  Properties  of  Carbon  Material  98 

66  Prediction  of  Bulk  Poisson’s  Ratio  of  Carbon  Foam  at  Various  Densities  98 

67  Prediction  of  Bulk  Tensile  Moduli  of  the  Carbon  Foam  Predicted  with 

Isotropic  and  Anisotropic  Properties  of  Foam  Ligaments  100 

68  Prediction  of  Bulk  Tensile  Moduli  of  the  Carbon  Foam  Predicted  with 

Isotropic  and  Anisotropic  Properties  of  Foam  Ligaments  100 

69  Variation  of  Bulk  Tensile  Moduli  of  the  Carbon  Foam  with  the  Increase 

of  Longitudinal  Modulus  of  Foam  Ligament  101 

70  Prediction  of  Bulk  Poisson’s  Ratio  of  the  Carbon  Foam  Predicted  with 

Isotropic  and  Anisotropic  Properties  of  Foam  Ligaments  101 

71  Prediction  of  Bulk  Poisson’s  Ratio  of  the  Carbon  Foam  Predicted  with 

Isotropic  and  Anisotropic  Properties  of  Foam  Ligaments  102 

72  SEM  Image  of  the  Mode  11  Crack  Tip  Region  at  Two  Consecutive 

Load  Levels  106 

73  SEM  Image  of  the  Mode  11  Crack  Tip  Region  in  the  Resin-Rich  Region 

of  the  Sandwich  Specimen  107 

74  Eully  Extended  Interlaminar  Crack  at  (a)  186- lb  Loading  and 

(b)  226- lb  Loading  108 

75  Compression  Eiber  Eailure  on  Inner  Surface  of  the  Bottom  Eace  Sheet 

at  226- lb  Load  109 


X 


LIST  OF  FIGURES  (Continued) 


Figure  Page 


76  Variation  of  the  Effective  Shear  Modulus  |i*  as  a  Function  of  a 

Concentration  of  the  Inclusions  c  116 

77  The  Functions  and  their  First  Derivatives  versus  Argument  x\:  f2(xi) 

(solid  \m€),fi(xi)  (dot-dashed  line),//  ’(xi)  (dotted  line)  118 

78  The  First  Few  Iterations  of  Statistical  Averages  «3ii>i(xi)\qxs\x^xi 

Estimated  for  the  Functions  f2(xi)  and  Nonstop  g(r)  hy  the  MEFM  119 

79  Variation  of  the  Effective  Modulus  L* 2222  as  a  Function  of  a  Coordinate 
X2  Estimated  by:  MEFM  and  Nonstep  Function  g(r)  (dashed  line), 

MEFM  and  Step  Function  g(r)  (dot-dashed  line),  Mori- Tanaka  and 
Nonstep  g{r)  (dotted  line),  Mori-Tanaka  and  Step  Function  g(r)  (dotted 

line)  126 

80  Variation  of  Statistical  Average  Stresses  Inside  the  Inclusions 
«322>i(x2)  (solid  line)  and  «3i2>i(x2)  (dot-dashed  line)  for  the  External 
Loading  at  the  Infinity  c°°ij  =  6i26j2  and  o°°ij  =  6ii6j2,  respectively,  and  a 

Step  Function  g  126 

81  The  Radial  Distribution  Functions  g(r)  versus  Relative  Radius  r/a 
Estimated  by  the  HCM  at  c=0.55  (solid  line)  and  c=0.5  (dotted  line), 

and  by  HCSM  at  c=0.55  (dot-dashed  line)  and  c=0.5  (dashed  line)  133 

82  Histogram  of  Fractions  p(nt)  of  Testing  Window  Containing  rit 
Inclusions  Generated  by  the  HCM  (solid  curve)  and  by  the  HCSM 

(dotted  curve)  134 

83  Length  of  Monte  Carlo  Simulation  versus  an  Area  Fraction  c 

Appropriate  to  the  HCM  (solid  curve)  and  to  the  HCSM  (dotted  curve)  135 

84  The  Radial  Distribution  Functions  g(r)  versus  Relative  Radius  r/a 
Estimated  by  the  IPSM  at  c=0.65,  N=3100:  150  Shaking  (solid  curve), 

100  Shaking  (dashed  curve),  30  Shaking  (dot-dashed  curve),  10  Shaking 
(dotted  curve)  136 

85  The  Radial  Distribution  Functions  g(r)  versus  Relative  Radius  r/a:  Single 
Sample  546-0181  (dashed  curve).  Averaged  Over  10  Samples  of 
546-0161  (dotted  curve),  at  c=0.65.  Averaged  Over  10  Samples  of 
546-0181  (dot-dashed  curve).  Averaged  Over  Eight  Materials  Produced 

by  the  Different  Technological  Regime  (solid  curve)  137 


XI 


LIST  OF  FIGURES  (Concluded) 


Figure  Page 

86  Radial  Distribution  Functions  g( r)  versus  Relative  Radius  r/a  Estimated 
by  the  Modified  CRM  at  c=0.60  (dashed  curve),  c=0.65  (dot-dashed 

curve),  c=0.70  (dotted  curve),  c=0.75  (solid  curve)  138 

87  Radial  Distribution  Functbns  g( r)  versus  Relative  Radius  r/a  Estimated 
by  the  Modified  CRM  at  A^=81 1  and  c=0.90  (solid  curve),  c=0.85 

(dotted  curve),  c=0.80  (dot-dashed  curve),  c=0.75  (dashed  curve)  138 

88  Radial  Distribution  Functions  g( r)  versus  Relative  Radius  r/a  Estimated 
by  the  Numerical  Simulation  (solid  curve),  from  Experimental  Data 
(dotted  curve),  by  the  Analytical  Approximation  [54]  (dot- dashed 

curve)  139 

89  Variation  of  the  Effective  Shear  Modulus  p.*  as  a  Function  of  a 

Concentration  of  the  Inclusions  c  140 

90  Variation  of  the  Relative  Effective  Shear  Modulus  as  a 

Function  of  a  Concentration  of  the  Inclusions  c  Estimated  by  the 
MEFM  with  g( r)  [54]  (dot-dashed  line),  by  the  MEFM  with  the  RDF 
Simulated  by  the  Modified  CRM  (solid  curve),  by  the  MEFM  with 
Well  Steered  g( r)  (dashed  curve),  and  by  the  MEFM  with 

Experimentally  Estimated  RDF  o  141 


xii 


LIST  OF  TABLES 


Table  Page 

1  Elastic  Properties  of  Constituents  6 

2  Applied  Stress  at  Debonding  and  Maximum  Stress  Coneentration 

Factors  as  a  Function  of  Off-Axis  Angle  10 

3  Material  Properties  Used  in  Analysis  19 

4  Stresses  in  Undamaged  IM7/5250-4  Composite,  AT=-205°C,  ac=0  20 

5  Stresses  in  Undamaged  IM7/5250-4  Composite,  AT=0,  Oc=l  20 

6  Comparison  of  Experiments  with  Predietions  Simulating  Various 

Failure  Modes  29 

7  Material  Properties  of  AS4/3501 -6  59 

8  Predicted  Effective  Young’s  Moduli  and  Poisson’s  Ratio  of  an 

Isotropic  Foam  94 


FOREWORD 


This  report  was  prepared  by  the  University  of  Dayton  Research  Institute  under  Air  Force 
Contract  No.  F33615-00-D-5006,  Delivery  Order  No.  0002.  The  work  was  administered  under 
the  direction  of  the  Nonmetallic  Materials  Division,  Materials  and  Manufacturing  Directorate, 
Air  Force  Research  Laboratory,  Air  Force  Materiel  Command,  with  Karla  L.  Strong 
(AFRL/MLBCO)  as  Project  Engineer. 

This  report  was  submitted  in  March  2002  and  covers  work  conducted  from 
30  November  2000  through  29  November  2001. 


EXECUTIVE  SUMMARY 


The  cruciform  specimen  design  was  modified  to  characterize  the  fiber-  matrix  interface 
under  combined  transverse  and  shear  loading,  which  allowed  a  mixed- mode  failure  envelope  to 
be  constructed.  The  micromechanical  failure  modes  in  a  unidirectional  IM7/5250-4  composite 
under  tensile  loading  parallel  to  the  fibers  to  identify  initiation  and  growth  of  microdamage  were 
analyzed  using  the  axisymmetric  damage  model  (ADM)  of  a  concentric  cylinder.  The  effective 
thermoelastic  moduli  of  a  damaged  multiphase  woven  oxide-oxide  composite  was  computed  and 
compared  to  experiment.  It  was  shown  that  additional  damage  mechanisms,  such  as  fiber- matrix 
interfacial  debonding,  delamination  or  effective  matrix  degradation,  could  lead  to  lower 
estimates  of  the  Young's  modulus  in  the  thickness  direction.  The  fiber  pushout  test  was 
described  using  solutions  to  several  boundary  value  problems  using  two  numerical  methods, 
namely  ADM  and  finite  element  method  (FEM).  The  study  illustrated  the  complexity  of 
micromechanical  stress  fields  at  the  fiber-matrix  interface  and  proposed  some  guidelines  to 
examine  in  these  key  boundary  value  problems. 

A  B-spline  analysis  method  (B-SAM)  with  overlapping  B-spline  approximation  functions 
of  displacements  along  with  a  curvilinear  transformation  was  developed  and  successfully  applied 
to  a  three-dimensional  analysis  of  a  woven  textile  composite  including  incorporation  of 
cruciform  unit-cell  specimens.  The  optimum  cruciform  geometric  configuration  to  generate 
failure  in  the  gage  section  was  determined  parametrically  and  confirmed  by  Moire 
interferometry. 

Accurate  stress  fields  of  laminated  composites  with  slanted  free  edges  near  ply  interfaces 
were  calculated  for  various  slanted  angles  based  on  the  B-SAM.  The  behavior  of  the  power  of 
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singularity  was  highly  dependent  on  the  slanted  angles  of  the  free  edge  with  excellent  agreement 
with  the  two -term  singular  asymptotic  solution. 

Carbon  foam  was  modeled  with  three-dimensional  microstructures  to  develop  a  basic 
understanding  of  the  performance  of  open-cell  foam  materials.  The  microstructural  properties 
were  then  correlated  with  the  bulk  properties  through  this  model,  which  should  provide  a  basis 
for  establishing  a  process-property  relationship  and  optimizing  foam  properties.  The  predictions 
indicated  that  the  effective  bulk  moduli  of  the  carbon  foam  were  dependent  more  on  the 
transverse  modulus  of  the  foam  ligaments  than  on  the  longitudinal  one,  while  the  Poisson’s  ratio 
was  dependent  on  both  the  longitudinal  and  transverse  moduli. 

Mode  II  cracks  were  observed  in  situ  using  a  high  resolution  Philips  XL30  E-SEM. 
Interlaminar  crack  tip  features  were  examined  in  unidirectional  fiber-reinforced  composites  with 
and  without  the  resin-rich  areas  at  the  midsurface  of  the  laminate  and  face  sheet  failure  in 
sandwich  plates.  Eurther  work  is  required  in  order  to  investigate  the  Mode  II  crack  extension 
mechanism  in  resin- rich  regions. 

Nonlocal  models  of  stress  field  concentrations  were  developed  to  examine  the  effective 
elastic  properties  of  random  structures  in  composites  (i.e.,  fibers  within  bundles).  The  iteration 
method  and  the  Eourier  transform  method  were  both  examined,  each  with  its  own  series  of 
advantages  and  disadvantages.  The  obtained  relations  depend  on  the  values  associated  with  the 
mean  distance  between  inclusions,  and  do  not  depend  on  the  other  characteristic  size.  The 
iteration  method  was  found  to  be  superior  to  the  Eourier  transform  method  in  all  cases  except 
when  an  analytically  explicit  relation  is  required. 

The  models  above  were  applied  to  functionally  graded  materials,  which  are  statistically 
inhomogeneous.  The  use  of  small  clusters  of  homogeneous  materials  allowed  an  analysis  of 
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these  statistically  inhomogeneous  materials.  This  can  be  considered  similar  to  the  analysis  of 
imbedded  polymer  crystals  in  an  amorphous  mass  or  the  percolation  of  chains  in  a  polymer. 

As  part  of  the  above  modeling,  the  random  organization  of  fibers  within  bundles  in 
composites  was  examined.  The  particle  reorganization  induced  in  a  computer  simulation  by 
shaking  is  subjected  only  to  geometrical  constraints,  whereas  for  real  structures  the  packing  is  far 
more  complicated.  Our  simulation  technique  is  able  to  isolate  the  fundamental  geometrical 
constraints  from  other  physics- mechanical  and  chemical  effects,  and  the  results  provide  a 
valuable  benchmark  for  evaluating  sophisticated  packing  schemes  used  to  model  real  composite 
materials.  While  taking  a  single  point  probability  density  (volume  fraction)  into  account  can 
provide  a  rough  estimation  of  bounds  of  effective  properties,  more  informative  characteristics  of 
the  point  set  are  obtained  using  statistical  second-order  quantities  (such  as  two -point  probability 
density,  second-order  intensity  function,  radial  distribution  function,  and  nearest  neighbor 
distribution),  which  examine  the  association  of  a  point  relative  to  other  points. 
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1.  MIXED-MODE  FAILURE  CRITERIA  USING  CRUCIFORM  GEOMETRY 


The  single- fiber  cruciform  design  was  successful  in  eliminating  the  influence  of  free- 
edge  stresses,  which  were  present  in  transverse  testing  of  conventional  straight- sided  specimens. 
In  this  work,  the  cruciform  specimen  design  was  modified  [1]  to  characterize  the  fiber- matrix 
interface  under  a  combined  state  of  transverse  and  shear  stresses.  This  was  achieved  by  utilizing 
a  cruciform  specimen  in  which  the  wings  of  the  cruciform  specimen  were  not  perpendicular  to 
the  loading  direction.  Instead,  the  ratio  of  the  normal  to  shear  loading  was  governed  by  the 
amount  of  angle  the  fiber  (oriented  along  Z')  made  with  the  bading  direction  (X-axis),  as  seen  in 
Figure  1 .  Note  that  this  figure  is  a  schematic  of  the  cross  section  passing  through  the  midplane 
in  the  thickness  direction  of  the  off-axis  cruciform  test  geometry.  The  fiber  extends  across  the 
full  width  (21)  of  the  sample  and  is  centered  with  respect  to  specimen  thickness  (t).  The  angle  0 
between  the  fiber  axis  and  the  loading  direction  is  referred  to  as  the  off-axis  angle. 
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Figure  1.  Cruciform  Specimen  with  Fiber  Oriented  at  Angie  6  to  the  Loading  Direction 

In  this  study  five  off-axis  angles,  namely  15,  30,  45,  60  and  75  degrees,  were  investigated 
in  a  model  system  consisting  of  140-pm-diameter  silicon  carbide  (SCS-0)  fiber  in  Epon  828 
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epoxy.  The  epoxy  resin  (from  Shell  Chemical  Co.)  was  cured  with  a  polyetheramine  (Jeffamine 
D-230  from  Texaco,  Inc.)  for  14  days  at  ambient  temperature.  Curing  at  ambient  temperature' 
eliminated  thermal  residual  stresses,  which  were  otherwise  induced  from  the  mismatch  in 
coefficients  of  thermal  expansion  of  fiber  and  matrix,  and  were  therefore  not  considered  in  this 
work.  Model  single- fiber  off-axis  cruciform  samples  were  made  using  standard  fiberglass  epoxy 
end  tabs  with  the  following  specifications:  wing  height  (2h)  equals  5.08  mm,  loading  arm  width 
(2a)  equals  7.62  mm,  sample  width  (2/)  equals  38.1  mm,  grip-to-grip  distance  (2g)  equals  50.8 
mm,  and  fillet  radii  of  curvature,  Ri  and  I^,  equal  7.62  mm  and  2.54  mm,  respectively. 

Initiation  and  growth  of  interface  debonds  were  detected  optically  by  observation  of 
variations  in  the  intensity  of  light  reflected  from  the  surface  of  the  fiber  during  loading.  The 
reflected  light  technique  [3]  has  been  extensively  utilized  to  detect  debonding  in  textile- sized 
fibers,  such  as  15-|J,m-diameter  Nicalon  fiber,  and  7-|im-diameter  AS4  graphite  fiber.  During  the 
test,  illumination  was  set  at  an  intensity  such  that  before  loading,  the  fiber  surface  was  on  the 
verge  of  being  shiny.  This  lighting  produced  a  bright  area  along  the  approximate  centerline  of  the 
fiber  that  was  present  throughout  the  test.  The  reflective  centerline  area  is  not  to  be  confused 
with  the  debond.  For  this  work  we  concentrated  our  efforts  on  measuring  the  stress  level  at 
which  debonding  initiates  in  off-axis  cruciform  specimens.  The  magnification  was  adjusted  to  a 
level  great  enough  to  allow  the  debond  length  to  be  determined  accurately  while  being  low 
enough  to  maintain  the  maximum  field  of  view.  The  tradeoff  between  field  of  view  and 
magnification  was  necessary  for  the  best  opportunity  of  capturing  the  initiation  of  debonding, 
since  the  entire  gage  length  could  not  be  imaged  at  once. 


'  Chemical  shrinkage  in  room  temperature -cured  epoxies  could  be  significant  [2],  However,  shrinkage- induced 
strains  were  neglected  in  this  study  since  no  data  was  available  for  Epon  828  cured  with  polyetheramine.  Current 
effort  is  underway  to  estimate  these  quantities  and  will  be  reported  elsewhere. 
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Figure  2a  is  a  photomicrograph  of  a  cruciform  specimen  with  the  fiber  oriented  at  45 
degrees  with  respect  to  the  loading  direction.  Note  the  white  line  in  the  central  region  of  the 
cruciform  sample  in  the  gage  area.  This  picture  was  taken  with  the  sample  illuminated  by  a  light 
source,  and  the  white  line  that  we  are  seeing  in  this  image  is  the  light,  which  is  reflected,  by  the 
debonded  fiber- matrix  interface.  Figure  2b  is  a  magnified  view  of  the  central  region  of  the 
sample  and  clearly  shows  the  illuminated  white  area  corresponding  to  the  debonded  interface. 
These  micrographs  in  Figure  2  are  therefore  a  visual  confirmation  of  fiber- matrix  debonding  that 
took  place  in  a  cruciform  specimen,  and  the  length  of  the  reflective  (white)  region  is  an 
approximate  measurement  of  the  extent  of  debonding. 


Debonded 

interface 


N., 


Figure  2.  (a)  Visual  Observation  of  Interfacial  Debonding  in  a  45-degree  Off-Axis  Specimen; 

(b)  Magnified  View  of  Illuminated  White  Region  Shown  in  (a). 


The  single- fiber  off-axis  cruciform  specimen  was  analyzed  using  3-D  FEA  employing  the 
ABAQUS  code  [4].  The  SiC  reinforcement  and  epoxy  matrix  were  modeled  with  3-D  eight- 


node  brick  elements  with  the  elastic  properties  listed  in  Table  1. 


Table  1 

Elastic  Properties  of  Constituents 


Constituents 

E,  GPa 

V 

SCS-0 

400.0 

0.15 

Epoxy 

3.44 

0.35 
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The  fiber- matrix  interface  was  assumed  to  be  perfectly  bonded,  since  we  wanted  to 
evaluate  the  stress  distribution  in  the  sample  prior  to  damage.  Due  to  the  large  number  of 
degrees  of  freedom  that  were  anticipated  (because  of  3-D  modeling),  meshes  were  generated 
such  that  computational  efficiency  would  be  maintained  without  sacrificing  accuracy.  Finer 
subdivisions  were  employed  in  the  regions  where  the  stress  gradient  was  expected  to  be  high, 
such  as  along  the  fiber-  matrix  interface  and  in  the  fillet  regions.  The  finite  element  (FE)  model 
was  loaded  by  applying  tension  along  the  X-axis  (as  shown  in  Figure  1)  by  means  of  constant 
displacement  of  the  end  nodes  to  simulate  clamped-end  conditions. 

The  variations  of  radial  (Ox  x  )  and  shear  (Ox  z  )  stresses  at  the  fiber- matrix  interface  along 
the  fiber  length  are  shown  in  Figures  3  and  4,  respectively,  as  a  function  of  off-axis  angle,  6. 
These  stresses  were  evaluated  in  the  X'-Z'  plane  located  at  Y=0  (i.e.,  through  the  x- section 
passing  through  the  midplane  in  the  thickness  direction).  For  parametric  analysis,  the  smaller 
radius  of  curvature,  F^,  was  additionally  reduced  to  1.524  mm  for  a  30-degree  specimen  and  to 
0.318  mm  for  a  15-degree  specimen.  This  reduction  was  done  in  order  to  maintain  a  circular 
fillet,  with  a  150-degree  rotation  for  the  30-degree  off-axis  specimen  and  a  165-degree  rotation 
for  the  15-degree  off-axis  specimen. 

As  expected,  the  radial  stress  concentration  factor  was  the  largest  for  the  90-degree 
specimen,  and  the  magnitude  decreased  as  the  off-axis  angle  0  was  reduced.  On  the  other  hand, 
the  shear  stress  at  the  interface  was  maximum  for  a  45-degree  off-axis  specimen,  as  seen  in 
Figure  4.  With  smaller  values  of  the  off-axis  angle,  it  was  observed  that  the  peak  in  the  shear 
stress  distribution  became  broader  and  flattened  out.  Thus,  even  though  the  maximum  shear 
stress  was  reduced  as  the  off-axis  angle  was  lowered  below  45  degrees,  significant  values  of 
shear  stresses  now  occurred  over  larger  fiber  lengths  in  the  central  region  of  the  cruciform 
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Normalized  distance  along  fiber  length 

Figure  3.  Variation  of  interfaciai  Radiai  Stress  Concentration  Factor  with  Off-Axis  Angie. 


Normalized  distance  along  fiber  length 

Figure  4.  Variation  of  Interfaciai  Shear  Stress  Concentration  Factor  with  Off-Axis  Angie. 

specimen.  Also,  the  location  of  the  maximum  shear  stress  changed  with  the  cmciform  geometry 
and  did  not  necessarily  occur  at  the  same  location  where  the  interfaciai  radial  stress  was  a 
maximum.  Moreover,  because  of  unequal  values  of  radii  of  curvature  (Ri  =  7.62  mm  and 
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R2  =  2.54  mm),  the  shear  stress  was  no  longer  antisymmetric  (the  two  shear  stress  peaks  had 
different  width  and  amplitude)  about  the  loading  axis  and  was  not  zero  at  the  specimen  center. 

A  minimum  of  eight  samples  were  considered  for  each  one  of  the  three  off-axis  angles 
considered,  namely,  30,  45  and  60  degrees.  Table  2  lists  the  average  value  (over  all  specimens) 
of  externally  applied  stress  at  debond  initiation  measured  for  the  specimens  tested  along  with  the 
maximum  radial  and  shear  stress  concentration  factors  evaluated  at  the  fiber-  matrix  interface  as  a 
function  of  off-axis  angle.  As  a  first  approximation,  the  maximum  stress  criteria  [5]  could  be 
employed  to  predict  failure  initiation  in  the  off-axis  cruciform  specimens,  since  we  had 
demonstrated  that  debonding  initiated  in  the  interior  region  of  the  specimen,  which  was  free  of 
initial  stress  singularities.  Both  the  initiation  and  growth  criteria,  however,  could  depend  on 
interaction  between  interfacial  shear  and  normal  stresses  and  their  respective  modes  of  energy 
release  rates.  Instead  of  utilizing  both  the  maximum  radial  and  shear  stress  concentration  factors, 
we  also  considered  the  following  two  alternate  options:  (i)  use  the  maximum  radial  stress 
concentration  factor  and  the  shear  stress  concentration  at  the  location  of  maximum  radial  stress, 
or  (ii)  use  the  maximum  shear  stress  concentration  factor  and  the  radial  stress  concentration  at 
the  location  of  maximum  shear  stress.  In  either  of  these  two  cases,  we  did  not  obtain  significant 
differences  (less  than  five  percent)  in  the  stress  concentration  values  compared  to  the  present 
approach.  There  was  yet  one  more  option  possible  and  that  was  to  use  the  radial  and  shear  stress 
concentration  factors  at  the  location  of  debond  initiation.  However,  precise  measurement  of  the 
failure  initiation  site  was  not  made  for  the  specimens  tested.  Thus,  in  the  absence  of  the  exact 
location,  the  maximum  values  of  stress  concentration  factors  served  as  useful  approximations. 
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Table  2 

Applied  Stress  at  Debonding  and  Maximum  Stress  Concentration  Factors 
as  a  Function  of  Off-Axis  Angle 


Off-Axis 
Angle  (°) 

Applied  Stress 
at  Debonding 
(MPa) 

Standard 

Deviation 

(MPa) 

Maximum  Radiai 
Stress  Concentration 
Factor 

Maximum  Shear 
Stress  Concentration 
Factor 

30 

41.37 

6.49 

0.3825 

0.7413 

45 

32.3 

3.40 

0.7741 

0.8529 

60 

31.7 

5.03 

0.9938 

0.7402 

90 

33.4 

2.3 

1.267 

0 

The  data  for  the  90- degree  transverse  specimen  listed  in  Table  2  were  taken  from  our 
earlier  work  For  that  90- degree  geometry  (with  equal  values  of  fillet  radii  of  curvature), 

microscopic  examination  of  the  fracture  surface  had  revealed  [6]  a  smooth  fiber  surface  which 
was  indicative  of  failure  occurring  due  to  normal  stress  alone.  Further,  the  radial  stress  was 
maximum  in  the  loading  direction  at  the  specimen  center,  while  the  shear  stress  was  zero  at  that 
location.  Those  stress  concentration  values  were  therefore  entered  in  Table  2  for  the  90-degree 
specimen.  In  contrast,  none  of  the  off-axis  cruciform  samples  tested  in  this  study  (30-,  45-  and 
60-degree  specimens)  failed  along  the  fiber  length.  The  average  applied  stresses  at  debond 
initiation  for  different  off-axis  angles  were  multiplied  with  their  corresponding  stress 
concentration  factors  and  the  results  plotted  in  the  normal- shear  stress  space,  as  shown  in  Figure 
5. 

A  quadratic  failure  envelope  given  by 

+  Cj  =  C2  (1) 

fit  the  experimental  data  rather  reasonably  well.  The  constants  ci  and  C2  were  evaluated  using 
the  measured  value  of  interfacial  normal  strength  and  the  45-degree  off-axis  cruciform  data.  The 
constant  ci  was  estimated  as  1.5,  while  the  shear  strength  of  the  interface  was  extrapolated  as 
34.5  MPa.  The  empirical  curve  given  by  Equation  (1)  was  found  to  be  a  good  fit  to  the 
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Figure  5.  Failure  Envelope  for  Off-Axis  Cruciform  Specimens. 

remainder  of  the  off-axis  data  (i.e.,  30-  and  60-degree  measurements).  Note  that  for  this  part  of 
the  study  we  had  neglected  the  cure  shrinkage -induced  stresses  since  no  data  were  available  for 
Epon  828  cured  with  polyetheramine.  Current  effort  is  underway  to  estimate  these  quantities, 
and  when  properly  accounted  for,  it  would  result  in  a  horizontal  shifting  of  the  failure  envelope 
to  the  left.  Nevertheless,  the  proposed  off-axis  specimen  design  promises  to  be  an  elegant 
method  for  extrapolating  the  shear  strength  of  the  fiber-  matrix  interface,  although  an  independent 
assessment  of  the  shear  strength  would  be  a  valuable  test  of  the  proposed  failure  criteria.  Other 
test  methods,  such  as  torsional  loading  of  a  shear  analog  to  the  cruciform- shaped  specimen,  or 
torsional  loading  of  a  rectangular  composite  /77,  are  currently  being  explored  to  evaluate  the 
shear  strength  of  the  fiber- matrix  interface,  which  will  be  subsequently  compared  with  the 
estimated  value. 
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1.1  SUMMARY  AND  RECOMMENDATIONS 


The  results  of  the  present  investigation  were  highly  encouraging,  since  a  simple 
polynomial  fit  to  the  data  from  two  different  angles  was  seen  to  predict  the  measurements  at 
other  off-axis  angles.  The  scatter  in  the  data  was  relatively  small,  while  the  analytical  model 
remained  quick  and  efficient  to  run,  similar  to  the  90- degree  cruciform  geometry  developed 
earlier.  Failure  initiation  took  place  in  the  interior  region  of  the  specimen,  which  was  free  of 
initial  stress  singularities.  Thus,  the  evaluated  strength  values  were  not  influenced  by  the  free- 
edge  effects  which  seem  to  dominate  some  of  the  more  common  test  methods,  such  as  slice 
compression,  microindentation,  pushout  and  fragmentation,  used  to  characterize  the  shear 
strength  of  the  fiber-matrix  interface. 
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2.  DAMAGE  EVOLUTION  AND  EAILURE  MODELING  IN  UNIDIRECTIONAL 

GRAPHITE/BMI  COMPOSITE 

In  this  study  an  integrated  experimental  and  analytical  approach  was  initiated  to  examine 
and  predict  the  micromechanical  failure  modes  in  a  unidirectional  composite  when  subjected  to 
tensile  loading  parallel  to  the  fibers.  Understanding  and  modeling  of  failure  in  advanced 
composite  materials  involves  several  different  and  interactive  failure  modes  within  the 
constituents  of  each  ply  and  between  plies  of  a  composite  laminate.  This  in  turn  results  in  a  need 
for  accurate  modeling  of  the  stress  distribution  within  the  plies  and  its  constituents  at  the 
micromechanical  scale.  An  investigation  of  the  initiation  and  propagation  of  damage,  such  as  an 
interfacial  debond  or  a  matrix  crack,  must  be  based  on  appropriate  failure  criteria.  Clearly  both 
the  initiation  and  growth  criteria  may  depend  on  interaction  between  shear  and  normal  stresses 
and  their  respective  modes  of  energy  release.  Numerous  failure  scenarios  are  therefore  possible 
given  the  right  conditions,  i.e.,  particular  combinations  of  relative  properties  of  the  constituents, 
ultimate  strengths  and  critical  energy  release  rates.  Our  objective  was  to  base  failure  on 
fundamental  properties  of  the  constituent  materials,  such  as  fiber,  matrix,  and  interface.  This 
approach  will  provide  the  basis  for  developing  a  fundamental  understanding  of  the  failure 
mechanisms  that  will  lead  to  appropriate  failure  criteria  for  a  unidirectional  ply  with  a  minimum 
number  of  tests. 

In  this  work  we  have  provided  direct  and  indirect  evidence  of  the  existence  of  a  number 
of  basic  modes  of  initial  damage.  Unfortunately,  it  was  not  possible  to  observe  directly  the 
development  and  growth  of  damage,  especially  that  which  existed  in  the  interior  of  a  sample. 
Rather,  surface  observations  were  taken  intermittently  throughout  the  loading  history,  and  optical 
and  scanning  electron  microscopy  was  performed  on  the  failure  surface.  This  provided  a  fairly 
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comprehensive  picture  of  the  failure  events  that  had  taken  place  but  left  unanswered  questions 
regarding  the  sequence  of  these  events.  These  questions  were  then  addressed  by  the 
micromechanics  models  [8,9]  to  postulate  logical  failure  scenarios. 

The  material  system  was  carbon  fiber- reinforced,  toughened  bismaleimide,  IM7/5250-4. 
This  system  exhibits  characteristics  ideal  for  a  generic  study  of  this  nature;  it  is  well 
characterized  and  exhibits  consistent  properties.  A  unidirectional  panel  was  laid  up  from  the 
unitape  prepreg  and  was  cured  in  an  autoclave  in  accordance  with  the  manufacturer’s 
recommended  cure  cycle.  Nonperforated  Teflon  bleeder  plies  were  used  to  make  smooth 
finished  surfaces  of  the  specimen.  The  cured  composite  panel  was  postcured  for  six  hours  at 
230°C  in  an  oven.  Rectangular  coupons  were  cut  from  the  panels  using  a  diamond -impregnated 
saw  blade.  Specimens  were  12.5  mm  wide  with  a  gage  length  of  150  mm.  One  of  the  specimen 
surfaces  was  polished  using  2- micron  and  0.5- micron  polishing  powder  to  enhance  the 
microscopic  image  for  microdamage  detection.  Prior  to  testing,  the  polished  specimen  surface 
was  examined  under  a  microscope.  No  damage  was  observed  in  any  as-processed  specimen.  All 
specimens  were  then  stored  in  a  desiccated  cabinet  until  testing.  Fiberglass/epoxy  end-tabs  were 
bonded  to  the  specimens  prior  to  testing,  which  was  done  at  room  temperature. 

All  specimens  were  loaded  in  tension  at  a  strain  rate  of  0.02/min.  These  specimens  were 
subjected  to  initial  loading  ranging  from  60  percent  to  95  percent  of  the  ultimate  strength  and 
unloaded  for  microscopic  examination.  Acoustic  emission  (AE)  was  monitored  in  all  tests  to 
determine  onset  of  damage  and  subsequent  failure  processes  in  the  specimen  during  loading. 
Photomicrographs  with  400  to  600  magnification  were  taken  randomly  in  the  central  region  (~10 
mm  long)  of  the  polished  surface  to  record  damage  such  as  fiber  breakage,  matrix  microcracking 
and  fiber/matrix  debonding.  It  was  assumed  that  these  observations  were  representative  of 
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damage  accumulation  over  the  length  of  the  specimen,  since  it  was  physically  not  possible  to 
examine  the  entire  surface.  This  will  be  verified  subsequently  by  making  random  observations 
of  damage  accumulation  in  other  regions  at  a  later  time. 

Figures  6  and  7  show  a  typical  stress- strain  curve  and  AE  count,  respectively,  for  an 
IM7/5250-4  unidirectional  composite.  For  this  specimen  the  first  significant  AE  event  occurred 
at  9.34  kN,  probably  resulting  from  onset  of  damage.  However,  upon  unloading,  no  visible 
damage  was  detected  under  the  optical  microscope.  It  is  possible  that  microcracks  may  have 
closed  upon  unloading,  and  we  were  therefore  unable  to  observe  them.  Future  effort  will  be 
directed  to  examine  initiation  of  microdamage  under  in  situ  loading.  Further  increase  of  AE 
activity  with  applied  load  is  indicative  of  additional  damage  accumulation. 

Figures  8,  9,  11  and  12  show  photomicrographs  of  various  damage  mechanisms  observed 
during  loading.  These  photomicrographs  were  taken  from  a  specimen  that  was  loaded  to  a 
composite  stress  of  2380  MPa  (failure  stress  was  2586  MPa).  Numerous  fiber  breaks  were 
observed,  as  shown  in  Figures  8  and  9.  While  some  of  the  fiber  breaks  were  isolated,  as  shown 
in  Figure  8,  there  was  clearly  a  large  row  of  fiber  breaks  lying  along  a  single  plane,  as  shown  in 
Figure  9.  This  latter  feature  was  mostly  observed  close  to  final  failure  of  the  unidirectional 
specimen.  In  addition,  multiple  breaks  of  some  fibers  were  also  observed,  as  shown  in  Figure  8. 
This  suggests  that  stress  transfer  had  taken  place  in  some  regions  such  that  high  value  of  the 
redistributed  fiber  axial  stress  had  developed.  Figure  10  is  a  histogram  of  measured  distances 
between  multiple  fiber  breaks.  While  there  was  certainly  some  scatter  in  the  measured  values, 
the  average  distance  between  multiple  breaks  of  individual  fibers  was  approximately  75  microns 
(average  of  130  fiber  break  distances). 
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Figure  6.  Composite  Stress-Strain  Curve 
for  iM7/5250-4  Composite. 


Figure  7.  Acoustic  Emission  Count  as  a 
Function  of  Load  Level. 
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Figure  8.  Typical  Fiber  Breaks  Observed  During  Tensile  Loading  of  iM7/5250-4  Composite. 
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Figure  9.  Large  Row  of  Fiber  Breaks  in  a  Singie  Piane  Observed  Ciose  to  Faiiure  Load. 
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Figure  10.  Histogram  of  Measured  Distance  between  Multiple  Fiber  Breaks. 
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Figure  11.  Matrix  Cracks  Bridged  by  Fibers. 
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Figure  12.  Interfacial  Debonding  Accompanied  by  Fiber  Break. 


In  addition  to  fiber  breaks,  we  also  observed  matrix  cracks  normal  to  and  bridging  the 
fibers,  as  shown  in  Figure  11.  In  some  cases  it  was  seen  that  the  matrix  cracks  occurred  in  a 
single  plane,  which  was  different  from  the  plane  in  which  fibers  broke.  The  matrix  cracks  were 
rather  widely  distributed.  Finally,  we  also  observed  debonding  of  the  fiber- matrix  interface  in 
conjunction  with  fiber  break,  as  shown  in  Figure  12.  However,  not  all  broken  fibers  had  a 
debonded  interface  at  the  broken  fiber  end.  From  the  photomicrographs,  an  average  value  of 
debond  length  was  measured  as  24  microns. 

The  information  obtained  from  the  photomicrographs  just  described,  however  valuable, 
was  incomplete  since  we  did  not  know  the  applied  stress  levels  at  which  the  various  events  took 
place.  Table  3  lists  the  thermoelastic  properties  of  IM7  fiber  and  5250-4  matrix,  where  E,  v,  G 
and  a  represent  Young's  modulus,  Poisson's  ratio,  shear  modulus,  and  coefficient  of  thermal 
expansion,  respectively,  while  the  subscripts  T  and  A  on  the  material  properties  denote  the 
directions  transverse  and  parallel  to  the  fiber  axis,  respectively.  The  properties  of  the  matrix  were 
measured  in  the  laboratory.  The  fiber  axial  Young's  modulus  was  taken  from  the  Hercules  data 
sheet  /707,  while  the  remainder  of  the  properties  were  back- calculated  using  NDSANDS  [11  ] 
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Table  3 

Material  Properties  Used  in  Analysis 


Property 

Material 

IM7  Fiber 

5250-4  Matrix 

ET(GPa) 

19.5 

4.28 

EA(GPa) 

276 

4.28 

Vt 

0.7 

0.36 

Va 

0.28 

0.36 

Ga  (GPa) 

70 

1.57 

aiCIO'TC) 

3.8 

46.8 

aA(10'TC) 

-0.05 

46.8 

Ultimate  strength  (MPa) 

4224 

77 

Fracture  toughness  (MPa  m''"") 

- 

0.88 

in  conjunction  with  the  measured  unidirectional  ply  properties  (fiber  volume  fraction  equals  0.6 
and  fiber  diameter  is  4.5  mierons).  The  fiber  tensile  strength  was  also  measured,  as  shown  in  the 
histogram  in  Figure  13,  with  an  average  value  of  4.05  GPa.  Our  intention  was  to  apply  the 
mieromeehanies  models  [8,9,11],  using  the  input  data  in  Table  3,  to  deduee  what  we  could  about 
the  initiation  of  damage  and  to  postulate  reasonable  failure  properties,  wherever  neeessary,  to 
carry  out  the  analysis  of  subsequent  failure  sequences. 

To  determine  the  damage  initiation  event,  we  caleulated  the  stresses  in  fiber  and  matrix 
corresponding  to  the  composite  stress  at  which  the  first  significant  AE  activity  took  place.  Tables 
4  and  5  list  the  eomputed  stresses,  using  the  generalized  plain  strain  model  [11],  under  a 
thermoelastie  eooldown  and  uniaxial  tension,  respectively.  Of  the  five  speeimens  tested,  the  first 
occurrence  of  significant  AE  activity  took  place  at  an  average  composite  stress  of  812  MPa.  At 
this  composite  stress  value,  the  fiber  and  matrix  axial  stresses  were  ealculated  as  1303  MPa  and 
76  MPa,  respectively  (including  residual  stresses).  Note  that  the  major  contribution  to  the  matrix 
axial  stress  came  from  the  residual  stress  problem. 

It  was  now  possible  to  conceive  of  two  different  seenarios  of  damage  initiation.  In  the 
first  case,  comparing  the  matrix  stress  (76  MPa)  with  its  strength  value  (77  MPa),  it  was  very 
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Figure  13.  Measured  Distribution  of  IM7  Fiber  Tensiie  Strength. 


Tabie  4 

Stresses  in  Undamaged  iM7/5250-4  Composite, 
AT  =  -205°C,  oc  =  0 


Radius,  gm 

Oz,  MPa 

09,  MPa 

Or,  MPa 

Fiber 

0 

-36.5 

-13.0 

-13.0 

2.25 

-36.5 

-13.0 

-13.0 

Matrix 

2.25 

54.6 

52.0 

-13.0 

2.905 

54.6 

39.0 

0 

Tabie  5 

Stresses  in  Undamaged  iM7/5250-4  Composite, 
AT  =  0,  Oc  =  1 


Radius,  gm 

Gz,  MPa 

09,  MPa 

Gr,  MPa 

Fiber 

0 

1.65 

-0.000507 

-0.000507 

2.25 

1.65 

-0.000507 

-0.000507 

Matrix 

2.25 

.0261 

0.00203 

-0.000507 

2.905 

.0261 

0.00152 

0 

likely  that  matrix  microcracking  would  take  place  at  or  close  to  the  corresponding  composite 
stress  level.  The  full-cell  crack  model  of  Pagano  and  Brown  [12]  was  used  to  simulate  an 
annular  matrix  crack  initiating  at  the  midpoint  between  fibers.  As  the  crack  grew,  the  tip  of  the 
matrix  crack  advanced  toward  the  fiber- matrix  interface.  We  found  that  the  potential  energy 
release  rate  at  the  matrix  crack  tip  attained  a  peak  value  as  the  crack  tip  approached  the  interface. 
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However,  the  model  predicted  that  there  was  insufficient  driving  force  to  mobilize  the  crack, 
since  the  calculated  peak  value  (1.1  J/rr?)  was  much  less  than  the  critical  value,  Qc  (157.5  J/nf ). 
The  latter  value  was  determined  via  the  relation 


K 


(2) 


where  Kic  was  the  fracture  toughness  of  the  matrix  material.  Consequently,  additional  damage 
had  to  occur  before  the  matrix  crack  became  unstable. 

We  also  examined  the  fiber  and  matrix  stresses  at  the  interface  as  the  matrix  crack 
approached  the  interface.  For  these  calculations,  the  fiber  radius  was  2.25  |a,m,  while  the  outer 
matrix  radius  was  set  at  2.905  |xm  in  order  to  obtain  60  percent  fiber  volume  fraction.  The  tip  of 
the  matrix  crack  was  successively  placed  at  a  radius  of  2.7,  2.5875,  2.43,  2.3175  and  2.25 
microns,  to  simulate  an  advancing  matrix  crack.  It  was  found  that  the  radial  stress  at  the 
interface  remained  largely  compressive  away  from  the  matrix  crack  plane,  while  the  shear  stress 
value  exhibited  a  peak  in  its  behavior,  as  shown  in  Figure  14.  The  axial  stress  in  the  fiber  was 
also  found  to  increase  in  magnitude  in  a  region  near  the  matrix  crack  plane  as  the  matrix  crack 
approached  the  interface,  as  seen  in  Figure  15.  It  is  to  be  noted  that  the  large  ADM  stress  values 
at  the  crack  tip  were  merely  approximations  of  elastic  singularities.  Away  from  the  plane  of  the 
matrix  crack,  the  fiber  axial  stress  recovered  its  far- field  value.  Thus,  secondary  damages  in  the 
form  of  either  shear  debonding  of  the  interface  or  fiber  breakage  in  the  plane  of  the  matrix 
microcrack,  were  possible  with  matrix  microcracking  as  the  initiating  damage  event.  In  the 
second  scenario  we  compared  the  fiber  axial  stress  (1303  MPa)  at  the  first  occurrence  of  a 
significant  AE  event  to  its  average  ultimate  strength  (4224  MPa).  The  use  of  maximum  stress  as 
a  failure  criterion  implied  no  fiber  breakage  at  the  composite  stress  level  (812  MPa),  even  though 
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Distance  from  crack  plane  (micron)  Distance  from  crack  plane  (micron) 

Figure  14.  Shear  Stress  at  the  Interface.  Figure  15.  Fiber  Axial  Stress  at  the  Interface. 


there  was  considerable  scatter  (standard  deviation  equals  894  MPa)  in  the  measured  strength 
values.  However,  when  the  composite  stress  was  increased  to  1300  MPa,  the  fiber  axial  stress 
increased  to  21 10  MPa,  which  was  well  within  the  scatter  band  of  fiber  failure.  Thus,  it  was  very 
likely  that  fiber  failure  would  take  place  at  the  higher  composite  stress  value. 

Next,  we  assumed  that  fiber  cracking  was  a  damage  initiation  event  and  put  a  penny¬ 
shaped  crack  in  the  fiber  in  the  analytical  model  [8].  Similar  to  the  matrix  microcracking 
problem,  we  observed  the  variation  of  stresses  at  the  fiber- matrix  interface  as  the  penny- shaped 
fiber  crack  grew  from  a  radius  of  1.35  pm  to  2.25  pm.  As  shown  in  Figure  16,  it  was  found  that 
the  interfacial  radial  stress  near  the  fiber  break  plane  was  tensile  and  increased  in  magnitude  as 
the  crack  advanced  toward  the  fiber- matrix  interface.  The  other  significant  stress  component 
was  the  axial  stress  in  the  matrix,  which  was  also  tensile  and  increased  considerably  with  crack 
growth,  as  seen  in  Figure  17.  Thus,  secondary  damages  in  the  form  of  either  tensile  debonding 
of  the  interface  or  matrix  microcracking,  were  possible  with  fiber  cracking  as  the  initiating 
damage  event.  The  potential  energy  release  rate  at  the  crack  tip  was  also  calculated.  However, 
the  calculated  values  could  not  be  compared  to  the  critical  energy  release  rate  since  no  fracture 
toughness  information  was  available  for  IM7  fiber. 
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Distance  from  crack  plane  (micron)  Distance  from  crack  plane  (micron) 

Figure  16.  Radial  Stress  at  the  Interface.  Figure  17.  Matrix  Axial  Stress  at  the  Interface. 

2.1  SUMMARY  AND  RECOMMENDATIONS 

The  work  presented  in  this  study  is  part  of  an  ongoing  effort  to  support  the  failure 
analysis  of  a  unidirectional  laminate.  Damage  observations  were  made  on  IM7/5250-4 
composite  at  some  selected  stress  levels.  We  calculated  fiber  and  matrix  stresses  at  the  first  sign 
of  damage  initiation  and  examined  in  some  depth  two  possible  damage  modes.  Some  initial 
predictions  were  made  based  on  the  micromechanics  model  in  conjunction  with  the  limited 
experimental  data  that  were  currently  available.  The  present  analysis  seemed  to  favor  matrix 
microcracking  (over  fiber  breakage)  as  the  initiating  damage  event.  However,  this  could  not  be 
verified  because  no  visible  damage  was  observed  in  the  microscope  upon  unloading  from  that 
stress  level.  The  proposed  failure  scenarios  need  to  be  validated  by  making  in  situ  damage 
observations.  Also,  the  validity  of  the  secondary  damage  modes  needs  to  be  established  by 
making  additional  observations  at  other  intermediate  composite  stress  levels.  Subsequent 
analytical  work  will  also  examine  debonding  of  the  fiber- matrix  interface  in  conjunction  with 
fiber  breakage  and/or  matrix  microcracking. 
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3.  MICROMECHANICAL  RESPONSE  OE  A  MULTIPHASE  COMPOSITE 


The  objectives  of  this  effort  were  to  support  the  development  of  a  revolutionary  ultra¬ 
compact  combustor  in  the  areas  of  material  selection  and  analytical  modeling.  Oxide-oxide 
ceramic- matrix  composites  consisting  of  oxide  fibers  embedded  in  an  oxide  matrix  are  currently 
the  subject  of  increased  interest  due  to  their  inherent  resistance  to  oxidation  and  excellent  tensile 
and  fatigue  properties  at  room  and  elevated  temperatures.  In  this  work  we  considered  2-D 
reinforcement  of  Nextel  720  fibers  in  an  eight-harness  satin  weave  architecture  in  a  matrix 
consisting  of  alumina  particles,  voids,  inherent  cracks  due  to  the  processing  conditions,  and  a 
silica  bonding  agent.  Predicted  quantities  included  the  thermomechanical  effective  moduli, 
thermal  conductivity  tensor  and  phase  stresses  due  to  processing  conditions.  The  effective 
composite  properties,  in  turn,  will  serve  as  input  parameters  in  a  3-D  FEM  to  analyze  the  stress 
fields  in  the  combustor  under  operational  conditions.  A  feedback  loop  between  the  FEM  and  the 
micromechanical  model  can  now  be  established  to  define  the  evolution  of  the  service- induced 
damage. 

The  approach  followed  a  building  block  scenario  in  which  the  microstructure  was 
homogenized  at  various  levels  depending  on  the  particular  quantities  being  sought.  Eirstly,  the 
effective  moduli  of  the  multiphase  matrix  material  were  computed  by  use  of  a  three-phase 
version  of  the  Mori- Tanaka  scheme  [13-15],  These  results  were  compared  with  multilevel 
representations  of  particles,  voids,  and  binder  possible  with  the  NDSANDS  model  [11], 
Secondly,  the  effect  of  yarn  crimp  on  the  3-D  effective  moduli  of  a  composite  were  determined 
by  comparison  of  a  fabric  reinforcement  model  [16]  with  that  for  a  cross-ply  laminate  and  also 
with  the  straight- fiber  micromechanical  stiffness  model  NDSANDS  [11],  Since  the  as- processed 
composite  was  known  to  contain  a  fairly  well  defined  distribution  of  cracks  normal  to  the  fiber 
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directions,  we  next  determined  the  effective  moduli  of  a  damaged  ply,  using  the  known  crack 
spacing  in  the  Schoeppner-Pagano  model  [17]  and  the  moduli  from  Pagano’s  3-D  exact  laminate 
elasticity  theory  [18],  Finally,  the  two  cracked  layers  were  assembled  and  the  moduli  computed 
from  the  exact  laminate  theory.  It  was  these  results  which  were  then  compared  to  experiment 
and  used  to  assess  the  quality  of  the  various  assumptions  invoked  in  the  calculations,  and  to 
define  the  influence  of  the  damage  on  the  effective  moduli.  Some  of  those  steps  are 
demonstrated  below  as  significant  results. 

For  the  material  system  under  consideration,  the  homogenized  matrix  consisted  of  the 
stiff  alumina  particles  finely  dispersed  in  silica  with  a  significant  amount  of  porosity.  We 
therefore  assumed  that  both  the  alumina  particles  and  the  pores  were  spherical  in  shape.  Using 
the  handbook  material  properties  of  alumina  (E=310  GPa,  v=0.15,  a=  8.5  e-6/°C,  p,=10  W/m-K) 
and  silica  (E=72  GPa,  v=0.16,  a=0.5  e-6/’C,  p=1.75  W/m-K),  the  effective  matrix  properties 
were  evaluated  as 

Em=  69.8  GPa,  Vm=0.18,  am=7.44  e-6/°C,  |Jm=2.66  W/m-K 
using  particle  and  void  content  of  0.52  and  0.34,  respectively,  in  the  matrix.  Secondly,  the  effect 
of  yam  crimp  on  the  3-D  effective  moduli  of  a  composite  was  determined.  The  homogenized 
matrix  moduli  were  assumed  in  this  comparison.  The  fabric  itself  was  an  eight-harness  satin, 
which  could  also  be  approximated  by  a  cross-ply  laminate.  Assuming  that  no  initial  damage  was 
present  and  using  ^262  GPa,  Vf=0.25,  af=6.0  e-6/°C,  pf=5.9  W/m-K,  for  Nextel  720  fiber,  we 
evaluated  the  following  properties  for  the  undamaged  fabric, 

Ex  =  Ey=  136.3  GPa,  116  GPa,  Vxy=0.1859,Vxz  =  Vyz=  0.2130 
Gxy  =  49.42  GPa,  Gxz=  Gyz=  48.75  GPa,  ax=  ay=  6.502  e-6/°C,  a^=  6.719  e-6/°C 
Px=  Py=  3.935  W/m-K  and  |lz=  3.752  W/m-K, 
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where  x-y  was  the  fabric  plane  and  z  was  the  thickness  direction,  and  the  fiber  volume  fraction 
was  0.45.  These  computed  thermoelastic  constants  using  NDSANDS  [11]  were  in  very  good 
agreement  with  the  analytical  results  of  Naik  [16] ,  where  the  actual  fiber  architecture  shown  in 
Figure  18  was  considered  in  a  unit  cell,  thereby  supporting  our  approximation  of  the  fabric. 

E,,  =  Ey  =  71  -  87  GPa,  E,  =  47  -  53  GPa,  Vxy  =  0.07  -  0.1 1,  Vxz  =  Vyz  =  0. 12,  Gxy  =  17  GPa, 
ttx  =  tty  =  5.85  -  6.73  e-6/°C,  =  6.99  e-6/°C  |ix  =  2.15  W/m-K,  |iz  =  1.62  W/m-K, 

so  that  the  effective  property  predictions  were  much  too  high.  However,  on  close  examination, 
extensive  damage  in  the  form  of  transverse  matrix  cracking  was  observed  as  a  result  of  the 
shrinkage  which  occurred  during  pyrolysis  processing  [19],  as  seen  in  Figure  19.  These 
microcracks  were  distributed  throughout  the  composite  with  spacing  ranging  from  50  to  200  |J,m 
prior  to  machining  and  testing  and  were  simulated  using  the  large  radius  axisymmetric  damage 
model  [17[. 

We  next  determined  the  effective  moduli  of  a  damaged  ply  using  the  known  crack 
spacing  in  the  Schoeppner-Pagano  model  [17]  and  the  moduli  from  Pagano’s  3-D  exact  laminate 
elasticity  theory  [18].  Finally,  the  two  cracked  layers  were  assembled  and  the  composite  moduli 
computed  from  the  exact  laminate  theory.  This  procedure  was  repeated  over  the  entire  range  of 
crack  spacing  measurements.  The  results  of  these  calculations  are  shown  in  Figure  20.  This 
modeling  approach,  which  incorporated  the  presence  of  transverse  cracking,  was  seen  to  produce 
very  good  agreement  with  the  experimental  measurements  for  all  of  the  effective  thermoelastic 
moduli,  except  Since  the  transverse  cracks  were  aligned  in  the  thickness  directbn,  they  did 
not  influence  the  modulus  in  that  direction,  which  remained  unchanged  from  the  previously 
reported  undamaged  value.  Traditionally,  a  more  simplified  approach  to  modeling  discrete 
damage  is  the  property  degradation  model  in  which  the  stiffness  of  the  load- carrying  element  is 
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1  mm 


200  |a,m 


Figure  18.  Microstructure  of  Woven 
Nextei  720/AS. 


Figure  19.  Magnified  View  of  Microstructure 
Showing  Matrix  Cracks  Due  to 
Processing. 


Figure  20.  Effective  Thermoeiastic  Moduii  of  Woven  Nextei  720/AS  Composite  as  a  Function  of 
Transverse  Crack  Spacing. 


reduced  to  reflect  the  presence  of  damage.  In  the  present  case  we  degraded  the  elastic  modulus 
of  the  matrix  material.  A  parametric  study  revealed  that  a  matrix  elastic  modulus  of  20  GPa 
(reduced  from  undamaged  value  of  69.8  GPa)  resulted  in  a  reasonably  good  match  with 
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experiments  for  composite  Young's  moduli  in  both  thickness  and  in- plane  directions,  as  shown  in 
Table  6. 

Next,  to  calculate  more  accurate  phase  stresses  so  that  the  failure  characteristics  could  be 
predicted,  we  replaced  the  previous  layer  models  with  a  concentric  cylinder  model  [9]  in  which 
the  distinct  fiber  and  homogenized  matrix  phases  were  recognized,  along  with  a 
micromechanical  representation  of  the  damage.  The  latter  damage  mode  consisted  of  fiber- 
matrix  debonding  (over  half  the  fiber  surface)  coupled  with  matrix  cracks  emanating  from  the 
debond  tips  (and  aligned  in  the  thickness  direction).  This  model  could  now  be  checked  for 
consistency  with  the  damage  observations.  As  shown  in  Table  6,  the  prediction  of  the  in-plane 
Young's  modulus  from  the  micromechanics  model  [9]  agreed  well  with  the  experimental  data. 

As  opposed  to  straight  transverse  cracks  in  the  thickness  direction  in  the  laminate  model  [17], 
the  micromechanics  model,  which  considers  180- degree  fiber- matrix  debonding  in  conjunction 
with  radial  matrix  cracking,  now  leads  to  a  reduced  modulus  in  the  thickness  direction  (as 
compared  to  the  undamaged  value).  The  predicted  value  of  ^  (72  GPa)  was  still  higher  than  the 
average  value  of  50  GPa  measured  in  the  thickness  direction.  To  realize  the  smaller  values  of  Ez 
which  were  measured,  we  once  again  appeal  to  the  layer  model  [17]  but  this  time  considered 
de lamination  between  the  0-  and  90-degree  plies  in  addition  to  transverse  cracking.  It  was 
assumed  that  ply  delamination  simulated  fiber-matrix  debonding,  which  could  occur  in 
conjunction  with  transverse  cracking.  The  predictions  of  this  simulation  are  also  shown  in  Table 
6  for  a  delamination  length  of  42  |a,m  and  transverse  crack  spacing  of  48  |J,m.  These  predictions 
of  Young's  moduli  seemed  to  reinforce  the  hypothesis  that  an  additional  damage  mode,  such  as 
fiber-  matrix  debonding  or  delamination,  was  occurring  in  conjunction  with  transverse  cracking, 
which  in  turn  leads  to  a  lower  value  of  ^  and  thus  a  better  match  with  the  measurements. 
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Table  6 

Comparison  of  Experiments  with  Predictions  Simulating  Various  Faiiure  Modes 


Ex  —  Ey,  GPa 

Ez,  GPa 

Experimental  Measurements 

71-87 

47-53 

Laminate  model 

{Undamaged} 

136 

116 

Laminate  model 

(T ransverse  cracking,  25  to  200  |im  crack  spacing} 

80-89 

116 

Matrix  Degradation 

{Em  =  20GPa} 

86 

43 

Micromechanics  model 

(Transverse  cracking  +  interfacial  debonding} 

86 

72 

Laminate  model 

(Transverse  cracking  (48  |im  crack  spacing) 

+  delamination  (42  jam  length)} 

78 

55 

Lastly,  Lu  and  Hutchinson  [20]  have  shown  that  matrix  cracking  in  combination  with 
interfacial  debonding  has  the  potential  for  signifieantly  redueing  the  overall  longitudinal  thermal 
conductivity  of  a  unidirectional  fiber- reinforced  composite.  Recently,  Islam  and  Pramila  [21] 
used  FEM  to  eonsider  combined  effects  of  partial  debonding  and  matrix  cracking  on  the  effective 
transverse  thermal  conductivity  of  fiber-reinforced  composites.  They  concluded  that  the 
reduction  of  effective  thermal  conductivity  may  be  as  large  as  50  percent  when  the  interfacial 
conductance  is  redueed  by  two  deeades.  In  the  absence  of  a  suitable  model  for  evaluating 
thermal  conductivity  of  the  damaged  fabric,  our  analytical  predictions  for  the  undamaged  fabric 
represent  upper  bound  results  and  are  therefore  signifieantly  higher  than  the  measured  values  for 
the  damaged  fabric. 

In  summary,  we  have  demonstrated  that  discrete  damage  modeling  was  essential  for 
effeetive  property  estimation.  From  the  meehanics  viewpoint,  we  have  provided  a  set  of 
analytical  tools  and  experimental  protocol  to  enact  true  composite  material  design  for  high- 
temperature  applications.  Aside  from  the  oxide-oxide  elass  of  eomposites  (or  possibly  another 
form  of  ceramie- matrix  eomposite),  the  methodology  is  also  appropriate  for  application  to 
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structural  carbon- carbon  composites  which  have  numerous  functions  in  thermal  protection 
systems  and  other  space  vehicle  applications. 
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4.  SINGLE-FIBER  PUSHOUT  CONSIDERING  INTERFACIAL  FRICTION  AND 

ADHESION 


We  completed  work  on  this  topic  and  published  a  journal  paper  [22]  on  this  subject.  The 
paper  presented  an  effort  to  understand  the  complexity  of  stress  distributions  and  deformation 
fields  at  the  fiber- matrix  interface  in  a  typical  single- fiber  pushout  specimen  subjected  to  a 
variety  of  interface  conditions.  Two  numerical  methods,  a  stress-based  Reissner's  variational 
solution  and  the  displacement-based  FEA,  were  used  to  compute  the  numerical  solutions.  The 
fiber-  matrix  interface  conditions  considered  included  zones  with  perfect  bonding,  frictional 
contact  (stick/slip)  and  open  crack  faces.  The  location,  size,  and  progression  of  the  interface 
zones  were  selected  to  mimic  realistic  (from  experimental  findings)  and  other  typical  interfacial 
damage  growth  during  the  pushout  test.  Load  cases  examined  include  the  residual  stress  due  to 
chemical  shrinkage  during  specimen  processing  and  the  applied  pushout  loading.  The  discussion 
presented  in  the  paper  focused  on  two  aspects  of  fiber  pushout  testing.  The  first  aspect  was 
extracting  information  pertaining  to  the  material  behavior  from  the  load-displacement  curves. 

The  second  aspect  was  deducing  the  nature  of  the  interface  from  the  test  results  with  the 
assistance  of  the  numerical  simulations. 

It  was  shown  that  the  stress  distributions  and  the  progression  of  interface  zones  largely 
depend  upon  the  initial  state  of  the  interface.  Comparisons  were  made  illustrating  the  differences 
in  the  stress  fields  due  to  adhesive  (no  relative  axial  displacement  throughout  the  loading  history) 
and  stick  (no  relative  axial  displacement  during  the  last  load  increment  but  may  have  nonzero 
total  relative  axial  displacement)  interface  conditions.  These  differences  could  become 
extremely  important  in  the  interpretation  of  analytical  results. 

It  was  further  pointed  out  that  although  knowledge  of  the  initial  damage  and  state  of  the 
interface  in  the  specimen  was  required,  the  presence  or  absence  of  adhesion  at  the  interface  could 
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be  inferred  and  characterized  based  on  numerical  simulations.  An  indicator  for  the  presence  of 
adhesion  at  the  interface  was  the  maximum  load  carried  by  the  specimen.  A  three- step 
methodology  was  presented  which,  when  used  in  concert  with  the  experimental  results,  could  be 
used  for  assessing  the  nature  of  the  interface  and  the  extent  of  adhesive  bonding.  However,  one 
could  not  define  the  detailed  zones  if  adhesion  was  present  unless  one  determined  the  debond 
length  throughout  the  test  [perhaps  by  the  nondestructive  evaluation  (NDE)  method]  or  had  a 
valid  debond  initiation  criterion  Some  of  the  suggested  simulations  extended  beyond  the  normal 
pushout  testing  protocol. 

Finally,  given  the  properties  and  the  initial  state  of  the  interface,  the  stress  distributions  at 
the  fiber- matrix  interface  were  shown  to  be  affected  by  the  size  and  order  of  the  interface  zones 
and  the  singularities  at  the  zone  transitions.  The  boundary  conditions  at  the  interface  changed 
continuously  with  pushout  loading  leading  to  a  continuous  change  in  order  and  size  of  interfacial 
zones.  A  series  of  stress  distributions  were  presented  under  several  interface  conditions  to 
illustrate  the  complexity  of  stress  distributions  during  interfacial  zone  growth  in  the  single- fiber 
pushout  testing.  These  stress  distributions  highlighted  that  neither  simplified  shear- lag  analyses 
nor  local  asymptotic  analyses  could  capture  the  stress  fields  that  drove  the  interfacial  damage 
development  in  this  specimen.  Accurate  global  analyses  were  necessary  to  capture  the  various 
mechanisms  in  progress  during  the  pushout  testing. 
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5.  PARAMETRIC  DESIGN,  ANALYSIS,  AND  TESTING  OF  A  MODEL  WOVEN 
COMPOSITE  WITH  A  CRUCIFORM  GEOMETRIC  CONFIGURATION 

To  achieve  the  optimum  structural  properties  of  state-of-the-art  fabric  reinforcements  of 
woven  composites,  there  is  a  need  to  develop  a  basic  understanding  of  deformation  and  damage 
mechanisms.  The  failure  in  textile  composites  in  general  occurs  in  the  vicinity  of  two 
perpendicular  yarns  crimping  over  and  under  each  other.  In  situ  experimental  observation  [23] 
of  damage  initiation  in  textile  composites  reveals  that  the  damage  initiates  in  the  form  of 
interface  cracks  in  the  vicinity  of  yarn  crimping,  which  is  strongly  influenced  by  the  interlaminar 
stresses  at  the  interface  region.  Thus  an  accurate  prediction  of  the  interlaminar  stresses  at  the 
interface  region  is  needed  to  reliably  analyze  damage  and  failure  in  woven  composites. 

The  B-SAM  is  a  previously  developed  computer  program  for  the  3-D  analysis  of 
anisotropic  laminated  composite  structures,  including  bonded  joints  and  bolted  joints.  The  B- 
SAM  has  the  advantage  of  an  efficient  modeling  and  solution  scheme  through  the  use  of 
overlapping  B- spline- based  shape  functions  [24,25].  The  overlapping  nature  of  the 
approximation  functions  keeps  a  minimal  bandwidth  profile  in  the  stiffness  matrix  of  the  system 
of  equations.  While  the  conventional  FE  (p-element)  has  Cf  continuity,  m-th  order  spline 
functions  have  continuity,  which  yields  continuous  derivatives  up  to  the  (m-7)-th  order. 
Therefore,  stress  and  strain  fields  as  well  as  the  displacements  are  evaluated  continuously  within 
the  homogeneous  domain  (subregions)  with  the  spline  method.  The  homogeneous  domain 
discretized  with  n  subdivisions  in  one  direction  needs  (n+m)  overlapping  m-th  order  spline 
functions.  Figure  21  shows  an  example  of  cubic  splines  overlapped  on  seven  subdivisions.  The 
spline  functions  are  evaluated  with  a  recursive  formula  developed  earlier  [26], 
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1 


7  Subdivisions 


0  0.143  0.286  0.429  0.571  0.714  0.857  1 

Figure  21.  Overlapping  Cubic  B-Spline  Functions  with  Seven  Subdivisions. 

The  B-SAM  was  modified  to  accommodate  the  curved  geometry  of  woven  composites 
with  a  curvilinear  transformation.  Two  tensor  transformations  were  needed  for  setting  up  a 
global  constitutive  relation  between  the  global  stresses  and  strains:  one  for  the  fiber  orientation 
along  the  yams  and  the  other  for  the  yam  crimping  with  respect  to  the  global  coordinate  system. 
The  first  one  was  for  the  coordinate  transformation  between  material  principal  axes  (on- axes) 
and  off- axes  following  yarn  crimping. 

The  B-SAM  was  applied  to  the  analysis  of  a  unit  cell  of  a  plain- weave  textile  composite 
to  examine  how  well  the  solution  satisfies  interlaminar  traction  continuity  conditions  at  the  yam 
interfaces.  The  B-SAM  was  also  used  to  analyze  a  model  woven  composite  with  1-D  yam 
crimping.  The  maximum  stress  failure  criterion  was  applied  to  the  stress  field  solution  to  predict 
the  initial  failure  strength  and  the  damage  mode  of  the  1-D  model  laminate.  As  observed  from 
experimental  and  numerical  analysis  of  the  model  laminate  with  a  straight  edge,  matrix  cracking 
was  initiated  at  the  free  edge  of  the  transverse  fill  yam  perpendicular  to  the  loading  direction 
[27],  Such  a  failure  mechanism,  caused  by  the  free  edge,  was  not  desirable  for  the  analysis  of 
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the  woven  composites  because  most  of  the  interest  lies  in  the  effect  of  yarn  crimping.  Therefore, 
to  eliminate  the  free-edge  effect,  cruciform  specimens  were  used  in  this  study. 

In  an  attempt  to  validate  analytical  predictions.  Moire  interferometry  was  employed  to 
obtain  an  accurate  full- field  distribution  of  in- plane  strain  on  the  face  of  the  cruciform  specimen. 

5.1  NUMERICAL  RESULTS  AND  PARAMETRIC  STUD  Y 
5.1.1  Unit  Cell  of  Plain- Weave  Composite 

The  B-SAM  was  applied  to  the  analysis  of  a  unit  cell  of  a  plain- weave  textile 
composite.  The  unit  cell  of  the  model  was  divided  into  several  homogeneous  subregions; 
each  subregion  was  occupied  by  a  characteristic  fabric  yarn  or  a  matrix,  as  shown  in  Figure  22. 
and  are  the  length  of  the  unit  cell  in  the  x-  (warp)  and  y-  (fill)  directions,  and  and 

half  of  the  thickness  of  the  warp  and  fill  yarns,  respectively.  The  yams  were  assumed  to  be 
transversely  isotropic,  and  the  matrix  was  assumed  to  be  isotropic.  Each  yarn  and  matrix 
subregion  of  the  unit  cell  was  discretized  into  several  subdivisions  in  the  longitudinal  and 
transverse  directions.  Cubic  splines  were  used  as  the  approximation  functions  of  the  unknown 
variables  in  the  x-,  y-  and  z-directions.  A  boundary  condition  applied  to  the  unit  cell  was  as 
follows: 

u(0,y,z)^0  ,  u(L^,y,z)^Ug 

v{x,0,  z)^  v{x,  Ly,  z)^  0  (3) 

w(x,  y,  0)  =  0 

Analysis  was  carried  out  to  calculate  the  stress  components  along  the  yarn 
boundaries.  The  stresses  calculated  in  the  global  coordinate  axes  were  transformed  into  the  local 
coordinate  axes  along  the  yarn  crimping.  Alternatively,  these  local  stress  components  along  the 
interfacial  surfaces  were  related  with  global  stress  components  by  the  slopes  of  the  interfacial 
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Figure  22.  Unit  Celi  of  a  Piain-Weave  Composite.  Numbers  in  circles  indicate  the  numbers  of 
subregions.  Subregions  1  and  4  represent  warp  yarns,  2  and  3  represent  fili  yarns, 
and  5  and  6  represent  matrix  materiai. 


surfaces  in  the  x-  and  y-directions  ( h^'^l  and  ),  where  k  denotes  the  subregion  number.  The 

normal  and  shear  stress  components  along  the  yarn  interface  boundary  were  thus  calculated  with 
the  following  relations: 


df  =af 


■  •  G 
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(k) 
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(k) 
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where  d.  ’  and  o/  '  are  the  stress  components  at  the  interface  of  the  k-th  subregion  in  local  and 

global  coordinate  axes,  respectively.  The  subscript  (k)  will  be  dropped  later  in  this  section  unless 
it  explicitly  refers  to  a  different  subregion. 

Figure  23(b-d)  shows  the  normal  (dj  =  d^ )  and  shear  (^4  =1"^  ,  dj  =z^ )  stress 

distributions  along  three  interfacial  lines  in  Figure  23(a).  The  stress  distributions  plotted  with 
triangles  and  circles  were  evaluated  from  the  upper  and  lower  subregions,  respectively.  The 
results  show  an  excellent  agreement  with  the  interfacial  normal  and  shear  stresses  from  the  lower 
and  the  upper  subregion,  which  indicates  the  present  method  is  a  reliable  way  of  calculating  the 
interlaminar  stresses.  Note  that  one  of  the  shear  stress  distributions.  Figure  23(c),  shows  a  large 
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(a)  Unit  cell 


(b)  Normal  stress  (<7  j, ) 


(c)  Shear  stress  (X  )  (d)  Shear  stress  (X^ ) 

Figure  23.  Interfacial  Normal  and  Shear  Stress  Distributions  of  a  Unit  Cell  of  a  Plain-Weave  Textile 
Composite. 


discrepancy  between  the  lower  and  upper  stress  values.  This  component  is  numerically  less 
significant  than  the  other  stress  components,  Figure  23(b,d),  when  comparing  the  absolute 
magnitude  of  the  stress  components.  Figure  23(d)  shows  a  shear  stress  jump  in  the  middle  of  the 
distribution  where  the  material  properties  are  discontinuous.  At  this  point  the  stress  has  weak 
singularity  due  to  the  material  discontinuity  and  thus  causes  a  high  stress  gradient.  The  reason 
for  the  stress  singularity  in  the  local  shear  stresses  is  that  the  stress  components  cannot  be 
determined  uniquely  because  of  the  zero  thickness  of  subregion  3  at  the  singular  point.  While 
two  subregions  (subregions  5  and  1)  were  considered  in  calculating  the  interfacial  stresses  at  the 
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left-hand  side  of  x  =  L^j  2 ,  three  subregions  (subregions  5,  3  and  1)  were  considered  in  the 
calculation  at  the  right-hand  side.  In  general,  it  was  impossible  to  satisfy  such  a  continuity 
condition  at  the  singular  point  with  polynomial  approximation  functions,  including  splines.  At 
this  singular  point  an  asymptotic  approach  should  be  used  to  determine  accurate  stress 
components. 

5.1.2  1-D  Model  Laminate  with  Cruciform  Geometry 

A  model  wo\en  composite  with  1-D  yam  crimping  was  analyzed  numerically  and 
experimentally.  The  B-SAM  yielded  accurate  and  continuous  stress  field  solutions  in  warp  and 
fill  yarn  subregions.  The  solution  also  satisfied  interlaminar  traction  continuity  conditions  at  the 
yam  interfaces.  The  maximum  stress  failure  criterion  was  applied  to  the  stress  field  solution  to 
predict  the  initial  failure  strength  and  the  damage  mode  of  the  1-D  model  laminate. 

As  observed  from  experimental  and  numerical  analysis  of  the  model  laminate 
with  a  straight  edge,  the  matrix  cracking  failure  was  initiated  at  the  free  edge  of  the  transverse  fill 
yarn  perpendicular  to  the  loading  direction.  Such  a  failure  mechanism,  caused  by  the  free  edge, 
was  not  desirable  for  the  analysis  of  the  woven  composites  because  most  of  the  interest  lies  in 
the  effect  of  yarn  crimping.  Therefore,  to  eliminate  the  free- edge  effect,  cruciform  specimens 
were  used  in  this  study.  Figure  24  shows  one  quarter  of  the  cmciform  geometry  that  is  cut  by 
symmetric  planes  and  meshed  for  numerical  analysis. 

The  general  configuration  of  the  cruciform  specimen  was  with  arms  loaded  in 
tension  and  arms  that  were  free  from  loading.  The  specimen  had  a  [90/0]s  stacking  sequence 
everywhere  except  along  the  centerline  of  the  unloaded  arms  where  the  stacking  sequence  was 
reversed  to  [0/90]s.  This  resulted  in  a  specimen  with  a  1-D  woven  configuration  across  the 
specimen  width,  extending  down  the  unloaded  arms.  The  purpose  of  the  cruciform  geometry 
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Figure  24.  Dimensions  and  Geometry  of  One  Quarter  of  a  Cruciform  Specimen. 

was  to  position  the  free  edge  of  the  woven  portion  of  the  composites  away  from  any  loading  (i.e. 
at  the  ends  of  the  unloaded  arms). 

Nontrivial  stress  components  were  taken  into  account  and  applied  by  the 
maximum- stress  failure  criteria  to  calculate  a  failure  index.  The  failure  index  ( k )  is  the  ratio  of 
the  stress  component  to  the  failure  strength  of  the  composite  material.  The  stress  concentration 
factors  were  calculated  by  taking  all  stress  components  at  the  bottom,  middle  and  top  surfaces  of 
each  yam  (lower  fill  [90°],  middle  warp  [0°]  and  upper  fill  [90°]). 

Figure  25  shows  concentration  factors  of  an  in- plane  longitudinal  stress  {G  j) 
along  the  middle  of  the  transverse  fill  yam  and  out-of- plane  shear  stress  (a  j )  along  the  middle 
of  the  longitudinal  warp  yam.  The  stress  concentration  factors  (kj  =  G jjY  ,  5 /S 13 )  were 

normalized  by  a  nominal  stress  concentration  factor  (k^  -G  ^jY),  where  G  ^  was  an  applied 
stress  on  the  fill  yam  away  from  the  yam  crimping  region,  and  Y  and  5;  5  were  transverse  and 
shear  strengths  of  the  fill  and  warp  yams,  respectively.  While  the  straight-edge  specimen  yielded 
the  highest  stress  concentration  at  the  free  edge,  the  cmciform  specimens  eliminated  the  free- 
edge  effect.  Furthermore,  the  cmciform  specimens  yielded  higher  shear  stress  concentration  at 
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1.5 


Figure  25.  Normalized  In-Plane  Longitudinal  and  Out-of-Plane  Shear  Stress  Concentration 
Factors  Along  the  Middle  of  Transverse  Fill  Yarn  in  1-D  Model  Laminate. 


the  longitudinal  warp  yarn  than  the  axial  stress  concentration  at  the  transverse  fill  yarn.  The  high 
shear  stress  concentration  may  cause  the  damage  at  the  yam  crimping  area  in  a  delaminating 
failure  mode,  which  may  increase  the  yam  crimping  effect  of  the  woven  composites. 

While  the  cmciform  specimens  eliminated  the  free- edge  problem  and  increased 
the  yarn  crimping  effect,  it  caused  a  stress  concentration  at  the  fillet  region  as  well.  A  parametric 
design  study  was  conducted  to  determine  the  optimum  cmciform  geometric  configuration  for  a 
1-D  model  woven  composite  to  identify  a  specimen  configuration  that  leads  to  failure  of  the 
composite  material  caused  by  the  textile  configuration  and  not  because  of  the  external  geometry 
of  the  specimen.  The  parametric  study  focused  on  producing  a  geometry  where  failure  was  more 
likely  to  occur  in  the  center  gage  section  of  the  specimen. 
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Figure  26  shows  the  ratio  of  stress  concentration  factors  in  the  fillet  region, 

(k)  fillet ,  to  those  of  the  woven  region  at  the  yam  crimp,  (k)^g^en-  The  stress  analysis  indicated 
that  the  transverse  tensile  stress  was  always  maximum  at  the  middle  surface  of  the  90°  ply  in  the 
fillet  region,  and  therefore  {k)fiiiet  =(5  JY  .  Meanwhile,  the  woven  stress  concentration  factor 

was  calculated  by  the  method  described  earlier,  and  therefore  {k)^e>ven  -  ^5  depending  on 

the  failure  mode.  Figure  26(a-c)  shows  that  the  fillet  stress  concentration  was  greater  than  the 
woven  stress  concentration. 
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Figure  26.  Ratio  of  Stress  Concentration  Factors  of  Woven  to  those  of  Fiiiet  with  Various 
Dimensions  of  the  Cruciform  Specimen  Geometry. 
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The  effect  of  the  fillet  stress  concentration  can  be  reduced  by  designing  the 
cruciform  specimens  by  varying  the  dimensions  of  the  cruciform  geometry.  The  fillet  stress 
concentration  reduced  rapidly  with  the  increase  of  the  fillet  radii,  as  shown  in  Figure  26(a,b). 
However,  as  the  fillet  radius  increased  further,  the  stress  at  the  yarn  crimp  also  became  less 
concentrated,  and  the  rate  of  the  reduction  of  the  woven  stress  concentration  factor  became 
greater  than  that  of  the  fillet  concentration  factor.  Therefore,  an  optimal  fillet  radius  can  be 
found  to  maximize  the  possibility  of  the  failure  due  to  the  yam  crimp.  The  same  effect  can  be 
achieved  by  decreasing  Lx2  or  increasing  Z.3,7,  as  shown  in  Figure  26(c,d).  However,  Lys  had  no 
influence  on  the  ratio  of  two  stress  concentrations. 

While  the  woven  stress  concentration  can  be  maximized  with  respect  to  the  fillet 
stress  concentration  with  optimal  dimensions  of  the  test  specimens,  it  can  be  further  affected  by 
varying  the  yam  waviness  ratio.  Figure  27  shows  a  result  of  the  parametric  study  by  varying  the 
yarn  waviness  ratio.  Different  yam  waviness  ratios  were  obtained  by  changing  in  Figure  24. 

At  low  waviness  ratio,  the  woven  composite  failed  by  transverse  matrix  cracking  either  in  the 
lower  or  upper  fill  yams.  As  the  waviness  ratio  increases,  out-of-plane  shear  failure  became 
dominant  at  the  middle  warp  yarn.  Furthermore,  at  approximately  a/?  =  0.087,  the  woven  stress 
concentration  exceeded  the  fillet  stress  concentration,  so  that  we  can  expect  the  initial  failure  in 
the  yam  crimp  rather  than  in  the  fillet. 

5.2  EXPERIMENT  AND  VALIDATION 

Test  specimens  were  prepared  by  embedding  the  model  yarn  crimping  in  a  composite 
plate.  After  drilling  four  holes  at  the  locations  of  the  fillet,  the  plate  was  cut  with  a  diamond  saw 
in  the  longitudinal  and  transverse  directions  at  points  tangential  to  the  hole  edges.  This  produced 
the  loaded  and  unloaded  arms  of  the  cruciform  specimen. 
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Figure  27.  Ratio  of  Stress  Concentration  Factors  of  Woven  to  those  of  Filiet  with  Various 
Yarn  Waviness  Ratio. 

The  geometry  of  the  test  specimen  was  set  based  on  the  above  parametric  study  as 
follows:  4;  =  =  L^2  ^3.8575,  =  10.0302,  L^4  ^15.9893,  =7.1402, 

Ly2  =  10.0302  ,  =  0 ,  R  =  9.525  (units  in  mm).  Half  of  the  thickness  of  the  warp  and  fill 

yarns  were  measured  by  micrographs  of  the  cross  section  as  =  0.128  mm  and  =  0.124 

mm,  respectively,  which  made  the  yam  waviness  ratio  be  ajX  =  0.063 .  Figure  28  shows  good 

agreement  between  the  measured  actual  geometry  and  the  modeled  geometry  expressed  with 
sinusoidal  functions. 

Moire  interferometry  [28]  was  employed  to  obtain  a  full- field  distribution  of 
displacements  and  strains  on  the  face  of  a  cmciform  specimen  of  optimum  geometry.  A  1,200- 
lines/mm  crossline  diffraction  grating  was  replicated  onto  the  face  of  the  cmciform  specimen. 
The  specimen  was  loaded  to  approximately  900  N,  and  displacement  fringe  patterns  were 
recorded.  Digital  phase  shifting  of  the  type  detailed  in  Lassahn,  et  al.  [29]  was  used  to  produce  a 
full- field  digital  array  of  displacement  values  over  the  central  region  of  the  specimen.  Numerical 
differentiation  was  used  to  calculate  the  in-plane  surface  strain  components. 
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Figure  28.  Comparison  of  Measured  and  Modeied  Cross-Section  Profiies  of  Modei  Woven 
Composite. 

Figure  29  shows  a  comparison  between  B-SAM  and  Moire  interferometry  results.  The 
B-SAM  results  were  scaled  to  match  the  far- field  axial  traction  applied  in  the  experiment.  The 
Moire  interferometry  results  have  been  smoothed  slightly.  Also,  because  of  slight  bending  of  the 
specimen  during  the  experiment,  the  Moire  results  have  been  averaged  in  a  symmetric 
arrangement  about  both  the  x-axis  and  the  y-axis.  Thus,  the  close-up  Moire  results  presented  in 
Figure  29  are  symmetric  about  the  x-  and  y-axes. 

Good  agreement  between  the  experimental  and  numerical  results  is  clearly  evident.  The 
other  two  components  of  experimentally  measured  in-plane  strains  yielded  similarly  good 
agreement  with  the  B-SAM  results.  All  of  the  in-plane  strain  components  sho^\ed  rapid 
variation  and  high  magnitude  in  the  fillet  region.  However,  only  the  axial  strain  showed 
significant  variation  in  the  region  of  yarn  crimping.  Nonaveraged  Moire  results  from  the  fillet 
region,  not  shown  in  this  paper,  also  displayed  a  close  match  with  the  B-SAM  predictions. 

5.3  SUMMARY  AND  RECOMMENDATIONS 

A  B-SAM  with  overlapping  B-spline  approximation  functions  of  displacements  along 
with  a  curvilinear  transformation  was  developed  and  applied  to  a  3-D  analysis  of  a  woven  textile 
composite.  Special  interest  was  given  to  the  reliable  calculation  of  interlaminar  stresses  at  the 
interfaces  between  yarn/yam  and  yarn/matrix.  A  unit-cell  of  plain- weave  composite  and  a  model 
woven  composite  with  1-D  yarn  crimping  were  analyzed  numerically  to  calculate  the  stresses 
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Figure  29.  Axial  Strain  Comparison  between  B-SAM  Predictions  and  Moire  Interferometry 
Measurements.  (Note  that  the  Moire  interferometry  results  were  averaged  in  a 
symmetric  manner  about  the  x-  and  y-axes.) 


and  to  predict  failure  strength  and  damage  mode.  Cruciform  specimens  were  used  to  eliminate 
free-edge  stress  concentration  of  the  model  woven  composite  by  positioning  the  free  edge  of  the 
woven  portion  of  the  composites  away  from  any  loading. 

A  parametric  design  study  was  conducted  to  determine  the  optimum  cruciform  geometric 
configuration  for  the  1-D  model  woven  composite  to  identify  a  specimen  configuration  that  led 
to  failure  of  the  composite  material  caused  by  the  textile  configuration  and  not  because  of  the 
external  geometry  of  the  specimen.  The  parametric  study  focused  on  producing  a  geometry 
where  failure  was  more  likely  to  occur  in  the  center  gage  section  of  the  specimen. 
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Moire  interferometry  was  employed  to  obtain  an  accurate  full- field  distribution  of 


in-plane  strains  on  the  face  of  the  cruciform  specimen.  These  data  showed  good  agreement  with 
the  numerical  results  obtained  with  B-SAM. 
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6.  ASYMPTOTIC  STRESS  ANALYSIS  OF  LAMINATED  COMPOSITES  WITH 

SLANTED  FREE  EDGE 


As  lightweight  and  high  stiffness/strength  material,  laminate  composites  have  been  used 
in  many  structural  applications,  such  as  airplanes,  ships,  automobiles,  sporting  goods,  etc.  An 
optimum  design  of  the  composite  structures  is  based  on  an  analysis  of  deformation  and  stress 
behavior  of  the  composites.  A  prevailing  method  in  the  analysis  of  laminate  composites  is  a  ply- 
level  homogenized  model,  which  considers  the  lamina  as  homogeneous  orthotropic  material. 

The  ply  level  model,  however,  results  in  singular  stress  behavior  in  the  vicinity  of  the  ply 
interface  and  laminate  edge.  The  singular  stress  is  caused  by  material  discontinuity  of  the  plies 
when  the  stress  tries  to  satisfy  traction  continuity  conditions  at  the  laminate  free  edge. 

Many  attempts  were  made  to  solve  the  straight  free- edge  problems  of  laminated 
composites.  A  mixed  approximation  of  displacement  and  stress  components  was  proposed  by 
Pagano  [30,31  ]  using  Reissner’s  mixed  variational  principle.  Without  the  precise  singular  stress 
terms,  the  mixed  method  demonstrated  a  highly  reliable  calculation  of  the  stress  fields  including 
interlaminar  stress  components.  Wang  and  Choi  [32]  constructed  an  infinite  series  elasticity 
solution  for  composite  wedge  near  the  ply  interface  and  the  wedge  edge  based  on  Lekhnitskii’s 
complex  variable  stress  function.  The  pioneering  work  by  Wang  and  Choi  [32]  was  extended  by 
Folias  [33]  and  Wang  and  Im  [34]  to  carry  out  an  asymptotic  solution  for  interlaminar  stresses  at 
the  interface  around  an  open- hole  edge  of  a  composite  plate.  They  showed  that  the  singular 
stress  solution  at  the  ply  interface  around  the  curvilinear  hole  edge  is  equivalent  to  that  for  the 
straight  edge,  provided  that  a  ply  orientation  with  respect  to  the  direction  tangential  to  the  hole 
edge  remains  the  same  as  the  ply  orientation  in  the  straight  edge,  with  respect  to  the  edge 
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direction.  Power  of  stress  singularity  was  found  to  be  dependent  upon  the  location  at  the  hole 
edge. 

larve  [35]  developed  a  B-SAM  for  3-D  stress  solution  of  laminated  composites  with  the 
open  hole.  A  numerically  calculated  stress  solution  near  the  hole  edge  by  this  method  was 
compared  with  one  given  by  the  singular  term  from  the  asymptotic  solution.  At  the  singularity, 
the  polynomial  spline  approximation  did  not  capture  directional  nonuniqueness  of  singular  stress 
functions.  However,  it  was  observed  that  the  singular  term  of  the  asymptotic  solution  with 
appropriate  coefficient  and  constant  additive  terms  matched  the  full- field  spline  solution  to  a 
distance  of  approximately  one  half-ply  thickness  from  the  singular  point.  The  large  area  of 
agreement  suggested  a  hybrid  method  by  superposing  the  singular  term  with  the  polynomial 
approximation  to  determine  stress  field  more  efficiently  and  reliably  without  a  tremendously  large 
amount  of  mesh  refinement  near  the  singular  free  edge. 

The  present  work  extended  the  asymptotic  analysis  to  the  problems  where  the  straight 
free  edge  is  no  longer  perpendicular  to  the  ply  interface.  The  slanted  free  edge  can  easily  be 
found  on  composite  joining,  such  as  in  countersunk  open  holes,  scarf- lap  bonded  joints,  etc.  The 
asymptotic  solution  was  to  be  compared  with  the  full- field  spline  approximation.  A  good 
understanding  of  the  singular  behavior  near  the  slanted  edge  will  lead  to  the  development  of  the 
hybrid  method  for  the  efficient  and  reliable  stress  analysis  of  slanted  laminated  composite 
structures. 

6.1  ASYMPTOTIC  ANALYSIS 

For  an  asymptotic  analysis,  three  coordinate  systems  were  introduced  for  laminated  plies 
with  slanted  free  edges  as  follows  (see  Figure  30): 
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Figure  30.  Global  and  Local  Coordinate  Systems  in  Asymptotic  Analysis. 


1.  xyz  ■  a  global  axis  in  a  Cartesian  coordinate  system. 

2.  rst :  a  local  axis  in  a  Cartesian  coordinate  system. 

3.  ri\|/^  :  a  local  axis  in  a  cylindrical  coordinate  system. 

Note  that  the  local  axes  were  variant  with  the  angle  (v|/  ).  An  assumption  was  made  that  the 
laminate  is  infinitely  long  in  the  global  y  -direction  (or  local  s  -direction),  so  that  the  analysis 

was  focused  on  a  2-D  cross  section  on  the  global  xz  -plane  (or  local  rt  -plane).  A  singular  point, 
which  is  the  center  of  each  coordinate  system,  was  located  at  the  intersection  of  the  free  edge  and 
the  interface  between  ( s  )-th  and  ( -I- 1  )-th  plies. 

In  the  local  coordinate  systems  (rir)/^  ),  r]  is  a  dimensionless  parameter  defining  the 
distance  from  the  free  edge,  and  V|/  is  an  angle  defining  the  direction  in  which  the  singular  point 

71  71 

is  approached.  In  this  coordinate  system,  0  <\|/  <  — -i-(|)  in  the  upper  ply  and  -¥  -  0 

the  lower  one,  where  (|)  is  the  slanted  angle  at  the  free  edge. 

The  global  ( xyz )  and  the  local  Cartesian  coordinate  systems  ( rst )  were  related  with  the 
following  local  transformation: 
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x|  mj  0  -/Tj  [r 

"  y  >  =  0  1  0  “  5  > ,  (5) 

z]  Jii  0  nil  _[t 

where  =cos\|/  and  /ij  =  sin\|/ .  Meanwhile,  the  global  ( xyz )  and  the  local  coordinate 
systems  (r|\|/^  )  at  the  cross  section  were  related  with  the  following  local  transformation: 

x=ricos\|/  ,  y  =  -^  ,  z  -  =r|  sin\|/ ,  (6) 


where  r|  >  0  and  (-7t/2  +([))  <\|/  <  (%/2+^) .  For  an  arbitrary  function  F(r| 
write 


J7 

dx  ^ 


dy  “  ac  ‘  az " 


(7) 


where 


AfF  =  costit 


dF  sin  \|/  dF 
ar|  r|  a\|/ 


.  ^  .  dF  cost!/  dF 

A„F  =  sm\it  — + - — 

ar|  r|  d\\f 


Equations  (7)  and  (8)  are  similar  to  the  ones  in  the  earlier  asymptotic  analysis  for  the 
laminated  composites  with  open  holes  [35].  Therefore,  we  can  utilize  a  similar  approach  to 
obtain  solutions  for  the  singular  analysis.  The  Navier  equations  of  equilibrium  for  a  given  ply 
can  be  written  as 


where  the  symmetric  square  (3  x  3)  matrices.  A,  B,  and  C,  are  expressed  as  follows: 

'a,  a.  0]“’  r  «  o  e.j+e,,]'''  [as  e«  oi"’ 

A=  e«  0  ,B=  0  2,6+845  .  C=  2.4  0  ,  (10) 

symm  ^55  symm  0  symm  ^33 
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where  are  off-axis  stiffness  coefficients  for  the  5-th  ply.  A  contracted  notation  is  used  for 
the  stiffness  and  stress  components  if  not  identified  otherwise,  so  that  indices  (/,  j  =  1,A  ,6 ) 
correspond  to  xx,  yy,  zz,  yz,  xz,xy ,  respectively. 

According  to  the  earlier  asymptotic  analysis  [35],  a  homogeneous  solution  of  the 
displacement  field  of  the  Navier  equations  in  Eq.  (9)  was  found  in  the  form  of 

,  i  =  1,2,3  ,  (11) 

k=l 

where  //*’  and  X  are  arbitrary  complex  constants  to  be  determined  from  the  boundary 
conditions.  Characteristic  values,  pf  * ,  are  eigenvalues  of  the  following  sixth-order  polynomial 
characteristic  equation, 

det| Ap^ -I- Bp  ^ -I- C  =  0  (12) 

and  dlf  are  eigenvectors  of  the  following  characteristic  matrix, 

dn 

(Ap^^ -I-Bpj^ -I-C)  ^  =  0.  (13) 

It  will  be  understood  that  coefficients  Pj.  and  vectors  are  constant  for  each  ply;  therefore,  the 
superscript  s  is  subsequently  omitted  unless  needed  for  clarity. 

For  orthotropic  materials  Eq.  (12)  yields  six  different  roots  such  that  p^  =  -p<-+3,  k  =  1,2,3 . 

For  isotropic  materials  there  are  only  two  different  roots,  ±  J-A  of  triple  multiplicity.  Since  the 
solution  in  Eq.  (11)  is  valid  only  if  all  characteristic  roots  are  distinct,  the  isotropic  material  is 
treated  as  a  slightly  orthotropic  material  by  perturbing  the  elastic  properties,  E  and  v  ,  with  a 
coefficient  of  anisotropy  (t  =10“’)  such  that 
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(14) 


E^={\-\-x)E  ,  E^=E^=E, 
^12=^13=^’  V23  =  (1  +  10t:)v  , 

GB  =  „,,f  G,,  =  G,3=(1+x)G„ 

2(1  +  V23) 

The  stress  components  in  the  global  ( xyz )  coordinates  are  expressed  as 


or  in  the  contracted  notation. 


(5  -  -■ 

IJ  ijkl  kl 


C.. 


ijkl 


du 
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V  dx,  ^ 


=  2  e 

J 


(15) 


(16) 


where  the  stress  components  a ,  correspond  to  ct  ^  ,a  ^  ,c>  ,a  ,c>  ,a  ^  ,  and  the  strain 

components  e,  correspond  to  ’^^yz  ’^^xy  *  =  hA  ,6,  respectively.  By 

substituting  Eq.  (11)  into  Eq.  (15),  one  can  obtain  the  global  stress  components  as  follows: 

6 

Gj  —  A,T|  (sin  \|t  +  cos\|r )  [Mit(2,i  2,6  ^112)'*’ Gis  ^113]  >  /  —  I,  2,  3,  6 

k=l 

6 

fk  (sin  ¥  +  Hi  costi/  )^-'  [2,5  d,,+  2,4  +  Hi  2,5  d,,  ]  ,  i  =  4, 5 

k=l 


The  power  X  and  the  coefficients  /j.  are  to  be  determined  from  the  boundary  conditions 


at  the  interface  between  the  two  plies  and  at  the  free  edges.  The  displacement  and  traction 
continuity  conditions  at  the  interface  between  the  5-th  and  ( 5  -I- 1  )-th  plies  yield 


~  K 

y  y  ~ 

Z 

^(s)  ^  (s+l) 

zz  zz  ’ 

,  0 

is)  _ 

yz 

(s+l) 


.  (s+1) 


at  x)/  =  0, 


(18) 


and  the  traction- free  boundary  conditions  at  the  slanted  free  edges  in  5-th  and  (  5-1-  l)-th  plies 
yield 
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(19) 


atX|/=-|+(^ 

The  stress  components,  {a  ,  in  the  local  coordinates  ( rst )  are  obtained  from  the  global 
stress  components,  {j  ,  as  follows: 

mf  0  nj  0  0 

0  1  0  0  0  0 

nl  0  ml  0  -  2m^n^  0 

0  0  0  mj  0  -nj 

—  m^n^  0  mj/Tj  0  mf -nf  0 

0  0  0  Hj  0  n\ 

Equations  (18)  and  (19)  form  a  homogeneous  system  of  12  linear  equations  with  12 
unknowns,  //*’  and  k  =1,A  ,6,  such  that 


M  f  =0 


where  M  is  a  coefficient  matrix  of  //*’  and  .  Nontrivial  solutions  of  Eq.  (21)  are  obtained 

by  requiring  the  determinant  of  the  coefficient  matrix  ( M )  to  be  equal  to  zero.  The  determinant 
of  the  system  of  equations  becomes  a  transcendental  equation,  which  yields  an  infinite  number  of 
roots  for  the  parameter  X  .  Note  that  only  roots  with  0  <  Re(A,)  <  1  provide  unbounded  stresses, 
which  dominate  the  solution  for  q  « 1 .  The  roots  for  the  determinant  were  found  numerically 
using  the  Secant  method.  Two  initial  guesses  of  X  were  used  to  obtain  the  determinants,  ^(Xq) 
and  g(X^) .  A  new  root  ( ^2 )  was  then  refined  by 


g(Xi)  (A,,  -)io) 

g(\)-gi\) 


(22) 
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After  setting  A.j  to  Xq  and  X2  to  A.j ,  the  procedure  repeated  until  X2  For  a  given  root  of 
X  ,  the  coefficients  and  ,  k  =  1,A  ,  6 ,  can  only  be  obtained  with  an  arbitrary 

multiplicative  factor.  Therefore,  one  of  the  coefficients  in  f  was  set  to  one. 

With  the  solution  of  Eq.  (21),  the  stress  components  in  Eq.  (15)  can  be  rewritten  as 

=L  ’  c,.(X.,¥)  >  j  =1,  A  ,  6,  (23) 

j=i 

where  ^  (Xj ,  V|/ )  are  determined  by  the  asymptotic  solution  and  normalized  so  that 
0^3  (Xj ,  -7U /2  +  (]))=  1 .  Even  though  the  stress  terms  in  Eq.  (23)  do  not  provide  the  full- field 
solution  of  the  3-D  problem,  the  singular  stress  terms  are  expected  to  correlate  with  the  full- field 
solution  for  r|  « 1 .  The  constants  (  Kj )  are  only  given  by  comparing  the  asymptotic  solution 

with  the  full- field  numerical  solution  for  q  « 1 .  As  stated  earlier,  the  singular  stresses  are 
obtained  only  for  roots  with  0  <  Re()i)  <  1 .  However,  when  the  constants  were  determined  by 
the  full- field  solution,  we  also  have  to  consider  the  roots  with  Re()i)  >  1  as  long  as  the  roots  are 
close  to  one,  because  of  weak  character  of  the  singularity  [35].  With  a  finite  number  of  the  roots 
(Xj  ,  j  -  I,  A  ,  n  ),  which  only  provide  meaningful  magnitudes  of  the  stress  components  for 
q  « 1 ,  Eq.  (23)  is  rewritten  as 

i  =  l,A,6,  (24) 

j=i 

where  the  additives  (F.)  are  also  given  by  comparing  the  asymptotic  solution  with  the  full- field 
numerical  solution. 
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6.2  FULL-FIELD  NUMERICAL  CALCULATION 


The  full- field  numerical  analysis  was  performed  with  a  laminated  composite  plate  of 
orthotropic  plies  having  a  slanted  edge,  as  shown  in  Figure  31.  The  slanted  angle  from  the 
global  z  -axis  was  denoted  as  (|) .  The  thickness  of  each  ply  was  denoted  as  h ,  and  other 
dimensions  in  x-  and  y  -directions,  when  measured  along  the  interface  between  the  plies,  were 
denoted  as  and  ,  respectively.  Total  thickness  of  the  laminated  plate  is  denoted  as  H  . 


The  plate  was  loaded  with  a  constant  displacement  (  m")  in  the  y  -direction  under  the  following 
boundary  conditions: 


Mj,(x,0,  z)  =  0,  Uy{x,Ly,z)  =  u‘' 

(L^,  y,  z)  =  (x,  y,  0)  =  (x,  y,  H)  =  0. 


(25) 


Zero  values  of  and  displacements  restrained  the  movement  of  the  planes  in  their 

perpendicular  directions,  which  simulated  the  plane  of  symmetry. 

The  full- field  3-D  solution  was  obtained  based  on  spline  approximation  of  displacement 
and  interlaminar  stress  components  in  curvilinear  coordinates.  The  stress  distributions  recovered 
near  the  free  edge  and  interface  between  plies  were  to  match  the  closed- form  asymptotic  solution 
results.  The  B-SAM  is  a  previously  developed  computer  program  for  the  3-D  analysis  of 
anisotropic  laminated  composite  structures,  including  bonded  joints  and  bolted  joints.  The  B- 
SAM  has  the  advantage  of  an  efficient  modeling  and  solution  scheme  through  the  use  of 
overlapping  B- spline- based  shape  functions  [36],  Furthermore,  the  spline  approximation 
eliminates  artificial  discontinuities  of  stress  and  strain  components  within  homogeneous 
domains. 

The  overlapping  nature  of  the  approximation  functions  keeps  a  minimal  bandwidth 
profile  in  the  stiffness  matrix  of  the  system  of  equations.  While  the  conventional  FE 
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(p  -element)  has  C°  continuity,  m-th  order  spline  functions  have  continuity,  which  yields 

continuous  derivatives  up  to  the  (m-l)-th  order.  Therefore,  stress  and  strain  fields  as  well  as  the 
displacements  are  evaluated  continuously  within  the  homogeneous  domain  (subregions)  with  the 
spline  method  for  m>2.  The  homogeneous  domain  discretized  with  n  subdivisions  in  one 
direction  needs  (n+m)  overlapping  m-th  order  spline  functions.  Figure  32  shows  an  example  of 
cubic  splines  overlapped  on  seven  subdivisions.  Construction  of  a  set  of  the  spline  functions  was 
achieved  with  a  recursive  formula  developed  earlier  /^557. 

The  B-SAM  analysis  needs  to  define  the  given  homogeneous  domain  in  X- space  with 
continuous  functions  in  ^  -space  defined  by  a  unit  cube,  so  that  any  variables  and  their 
derivatives  were  determined  from  the  functions  continuously  within  the  domain.  Figure  33 
shows  a  cross  section  with  slanted  edges  on  the  xz  -plane  mapped  to  the  ^  -domain. 

The  coordinates  in  the  X- domain  were  correlated  with  the  ones  in  the  ^  -domain  as 
follows: 


4  4 

a=\  a=\ 


(26) 
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7  Subdivisions 
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Figure  32.  Overlapping  Cubic  B-Spline  Functions  with  Seven  Subdivisions. 


Figure  33.  Mapping  of  Homogeneous  Domain  in  X-Space  to  a  Unit  Cube  in  ^  -Domain. 


where 

A',  =  (1-4.) (1-4,).  =4.(1-^,).  N,=iX.  (27) 

Note  that  Eqs.  (26)  and  (27)  are  valid  only  for  the  geometry  with  a  uniform  cross  section  along 
the  y  -axis.  The  slanted  geometry  in  different  planes  should  be  used  with  a  different  formula. 

The  variational  procedure  in  the  B-SAM  yielded  a  system  of  equations  for  spline 
approximation  coefficients  of  the  displacement  components.  Strains  and  stresses  were  then 
directly  calculated  with  the  derivatives  of  the  displacement  approximation  coefficients.  All 
stress  components  presented  in  this  section  were  normalized  with  an  average  stress  (o^ ).  The 
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average  stress  was  calculated  with  the  stress  compone  nts  in  the  loading  direction  on  the  cross 
section  perpendicular  to  the  loading  direction,  such  that 

A 

Quadratic  or  cubic  splines  used  for  displacement  approximation  were  continuously 
differentiable  at  all  points  in  the  ply.  Therefore,  all  strain  and  stress  components  were 
continuous  inside  each  ply.  Meanwhile,  the  interlaminar  strains  and  stresses  calculated  in 
adjacent  plies  may  be  discontinuous  at  their  common  interface.  However,  the  interlaminar  stress 
components  at  the  adjacent  plies  must  be  of  equal  value  with  each  other  to  satisfy  the  traction- 
continuity  condition  at  the  interface.  This  interlaminar  stress  continuity  condition  as  well  as  the 
zero  traction  condition  at  the  slanted  free  edge  will  be  investigated  by  comp  aring  with  the 
analytic  asymptotic  analysis  in  this  work. 

The  accuracy  of  the  numerical  solution  to  satisfy  the  conditions  described  above  was 
mesh  dependent.  Different  subdivisions  were  tried  to  determine  a  reasonably  refined  mesh  for 
the  given  geometry  and  the  material  properties.  Nonuniform  mesh  was  used  to  concentrate  the 
subdivisions  to  the  singular  point  at  the  intersection  of  the  free  edge  and  the  interface  between 
the  plies.  The  number  of  subdivisions  in  a  direction  in  the  domain  was  denoted  as  n .  The 
domain  was  subdivided  into  (n+1)  points  ( p, ,  /  =  0,  A  ,  n)  in  that  direction.  A  space  ratio  (q)  of 
the  nonuniform  mesh  was  given  such  that, 

q^PM-Pi  ^  /  =  (29) 

Pi-Pi-i 
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6.3  RESULTS 


A  symmetric  AS4/3501-6  laminated  composite  plate  with  a  lay-up  of  [45/-45]s  was 
considered  for  the  asymptotic  and  numerical  analysis.  The  dimensions  of  the  plates  are  =  10 
mm  and  t  =  0.125  mm,  and  the  material  properties  of  the  AS4/3501-6  ply  are  listed  in  Table  7. 


Table  7 

Material  Properties  of  AS4/3501  -6 


K  [GPa] 

[GPa] 

> 

II 

> 

=  G,,  [GPa] 

[GPa] 

139 

10.34 

0.3 

0.55 

5.52 

3.31 

The  asymptotic  analysis  was  first  performed  at  the  [45/-45]  interface  with  various  slanted 
angles  ((|) ).  The  characteristic  values  (  )  in  Eq.  (13)  were  ±  0.2827  i ,  ±  0.8581  i  and 

±1.209  i  for  both  45  and  -45  plies.  The  slanted  angles  varied  from  0  to  90  degrees  in  one- 
degree  increments.  As  stated  earlier,  the  roots  (.^)  were  found  in  the  range  of  0  <  Re()i)  <  1  and 
Re()i)  >  1  as  long  as  the  roots  are  close  to  one.  The  power  of  singularity  is  1  -  A, .  Distribution 
of  real  and  imaginary  parts  of  the  roots  with  respect  to  the  slanted  angle  is  plotted  in  Figure  34 
with  solid  and  dotted  lines,  respectively.  The  roots  were  pure  real  numbers  with  0°  <  (j)  <2°  and 
38°  <  (|)  <  90° ,  while  nonzero  imaginary  parts  of  the  roots  existed  with  3°  <  (|)  <  37°  .  When  the 
root  was  a  complex  number  with  the  non- zero  imaginary  part,  the  second  root  was  the  conjugate 
of  the  roots.  The  power  of  singularity  increased  significantly,  as  the  slanted  angle  increased  to 
be  greater  than  45  degrees. 

The  full- field  numerical  analysis  used  the  mesh  in  the  homogeneous  domains  with  the 
following  number  of  subdivisions  (m)  and  the  space  ratio  (q): 
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Figure  34.  Distribution  of  Reai  and  imaginary  Parts  of  the  Roots  with  Various  Slanted  Angles. 

1.  m  =  24,  ^  =  1.2  in  the  x-direction  for  both  plies. 

2.  m  =  6,q  =  0.5  in  the  y-direction  for  both  plies. 

3.  m  =  10,  ^  =  2.0  in  the  z- direction  for  lower  -45°  ply. 

4.  m  =  10,  ^  =  0.5  in  the  z- direction  for  upper  45°  ply. 

The  multiplicative  factors  ( Kj )  and  the  additives  ( F^)  in  Eq.  (24)  were  determined  by 

comparing  the  stress  components  from  the  asymptotic  solution  with  ones  from  the  full- field 
numerical  calculation.  In  later  figures,  the  stress  components  from  numerical  and  asymptotic 
solutions  will  be  plotted  with  solid  and  dotted  lines,  respectively.  The  stress  results  were 
compared  in  two  different  ways: 

1.  Stress  distributions  in  global  xyz  -axes  along  straight  lines  parallel  to  the  interface 
line  between  the  plies:  The  distance  between  two  lines  was  represented  with  a 
dimensionless  parameter,  z*  -  (z- z^'^)/h ,  where  z  and  z^*'  are  the  coordinates  of 
parallel  and  interface  lines,  respectively,  and  h  is  the  thickness  of  the  ply.  The 
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distance  was  chosen  at  z*  =  ±0,  z*  =  ±0.02  and  z*  =  ±0.04,  where  the  positive 
values  were  for  the  lines  in  the  upper  45  ply,  and  negative  for  the  lines  in  the  lower 
-45  ply. 

2.  Stress  distributions  in  local  rst  -axes  with  V|/  =(|)  along  contours  of  constant  radii 
with  the  center  at  the  singular  point:  The  radius  of  the  contour  was  denoted  as 
p  =  r|  .  The  radius  was  chosen  at  r]  =  0.02 ,  r|  =  0.04  and  r|  =  0.06 . 

Four  slanted  angles  were  selected  for  representing  distinguished  cases  as  appeared  in 
Figure  34,  such  that: 

1.  (|)  =  1° :  Roots  are  two  pure  real  numbers.  Both  roots  are  very  close  to  and  less  than 
one.  This  case  should  be  nearly  identical  to  the  problem  with  (|)  =0°  [35], 

2.  (|)  =  30° :  Roots  are  complex  conjugated  numbers.  The  magnitude  of  the  complex 
root  is  very  close  to  one. 

3.  (j)  =  45°  :  Roots  are  two  pure  real  numbers.  Both  roots  are  fairly  close  to  one,  and 
one  of  the  roots  was  greater  than  one. 

4.  (|)  =  60°  :  Roots  are  two  pure  real  numbers.  Both  roots  are  far  from  one,  and  one  of 
the  roots  was  greater  than  one. 

Note  that  the  roots  in  Figure  34  are  the  only  roots  whose  magnitudes  are  near  one,  whether  the 
roots  are  pure  real  numbers  or  not. 

4)  =1° 

Two  roots  found  by  the  asymptotic  analysis  were  =  0.9626  and  X2  =  0.9987 .  The 
second  root  was  practically  equal  to  one,  which  can  be  considered  as  a  nonsingular  root.  Figure 
35  shows  distributions  of  amplitudes  of  the  singular  stress  components  in  the  global  xyz  -axes  as 
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....... 

— 

(a)  =  0.9626  (b)  =  0.9987 


Figure  35.  Distribution  of  Magnitudes  of  Singular  Stress  Components  with  Two  Roots  for  (j)  =  1°. 


functions  of  \\f  for  the  two  roots.  Stress  components  a  and  o  are  nearly  constant  for  both 
roots.  The  interlaminar  stress  components  and  )  were  continuous  at  the  ply 

interface  (\|/  =  0 ),  which  satisfied  the  traction  continuity  boundary  condition  in  Eq.  (18).  With 
the  slightly  slanted  angle  of  the  free  edge,  all  stress  components  became  nonsymmetric  with 
respect  to  V|/  =  0 ,  which  were  otherwise  symmetric  if  (|)  =0°  [35], 

Figures  36  and  37  show  stress  distribution  by  asymptotic  and  numerical  analyses  in  the 
global  xyz  -axes  and  the  local  rst  -axes.  The  comparisons  in  both  figures  were  made  with  the 
same  multiplicative  coefficients,  Ki  =  -0.80  and  K2  =  0.08,  corresponding  to  the  first  and  the 
second  power  of  singularities.  As  stated  earlier,  the  multiplicative  coefficient  of  the  second  root 
was  much  less  than  that  of  the  first  root,  since  the  second  root,  which  was  nearly  one,  contributed 
minimally  to  the  asymptotic  stress  components.  The  additives  in  Figure  36  were  Fj  =1.15, 

F3  =  1.03  ,  F4  =  0.02 ,  Fj  =  -0.019  and  Fj  =  Fg  =  0 ,  and  in  Figure  37  were  F2  =1.15, 
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(c)  Interlaminar  normal  stress  (CJ 


(d)  Interlaminar  shear  stress  (CJ  ) 


(e)  Interlaminar  shear  stress  (G  ) 


(f)  In -plane  shear  stress  (G  ^ ) 


Figure  36.  Stress  Distribution  in  Giobai  xjz-Axes  by  Asymptotic  and  Numericai  Anaiyses  near 
the  Singuiar  Point  Aiong  Straight  Lines  Paraiiei  to  the  interface  between  Two  Piies. 
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Crt/co  Gtt^Go  Grr/c 


(c)  Interlaminar  normal  stress  (CJ^^) 


(d)  Interlaminar  shear  stress  (CJ  ) 


(e)  Interlaminar  shear  stress  (G 


(f)  In-plane  shear  stress  (G 


Figure  37.  Stress  Distribution  in  Locai  rst  -Axes  with  \|/  =  by  Asymptotic  and  Numericai 
Anaiyses  Aiong  Contours  of  Constant  Radii  with  the  Center  at  the  Singuiar  Point. 
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Fj  =  1.02 ,  F4  =  0.03  and  Fj  =  F5  =  Fg  =  0 .  Note  that  the  same  additives  were  used  for  both 
45°  and  -45°  plies. 

In  Figure  36,  with  the  range  of  x//i  >  0.01 ,  excellent  agreement  was  observed  between 
the  full- field  numerical  solution  and  two-term  asymptotic  one  for  all  stresses  with  one  set  of  the 
multiplicative  factors  and  the  additives.  One  exception  was  the  interlaminar  normal  stress  (<7^.^), 
which  showed  a  rapid  increase  in  the  asymptotic  solution  as  the  point  moved  away  from  the 
singular  point.  All  stress  components  grew  infinitely  on  approaching  the  singular  point  along  the 
interface  line  ( z*  =  ±0 )  in  the  asymptotic  solution,  but  yielded  finite  values  in  the  numerical 
solution.  Three  stress  components  ( and  d"^ )  from  both  analyses  were  nearly  equal  to 

zero  at  the  free  edge  (\|/  =  ±7t/2-i-(|) )  if  the  parallel  lines  were  not  along  the  interface.  Note  that 
these  stresses  in  the  asymptotic  solution  were  identically  zero  with  the  nonslanted  edge  problem 
(‘t>=0°). 

In  Figure  37,  three  stress  components  (d^^^ ,  d’^^  and  d”^,)  from  the  asymptotic  analyses 
were  identically  zero  at  the  free  edge  (\\f  =±7ij2+^  ),  which  satisfied  the  traction- free  boundary 
condition  in  Eq.  (19),  while  the  numerical  solution  yielded  nearly  zero  stresses.  Excellent 
agreement  was  also  observed  between  the  full- field  numerical  solution  and  two-term  asymptotic 
one  for  all  stresses  along  the  three  contours.  One  exception  was  again  the  interlaminar  normal 
stress  (cy„),  which  showed  more  discrepancy  as  the  radius  of  the  contour  increased.  This  bigger 
discrepancy  with  the  increase  of  the  radius  is  consistent  with  the  observation  made  in  Eigure  36, 
since  the  stress  increased  rapidly  in  the  asymptotic  solution  as  the  point  moved  away  from  the 
singular  point. 
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^  =30° 


No  real  root  was  found  near  one  with  this  slanted  angle.  Instead,  a  pair  of  the  conjugated 
complex  roots  was  only  found  by  the  asymptotic  analysis  with  =  0.9646  ±  0.0968  i .  The 

complex  roots  resulted  in  complex  numbers  of  the  magnitudes  of  stress  components  )  and  of 
the  multiplicative  coefficients  ( K. ).  Note  that  the  conjugated  roots  yielded  the  same  distribution 

of  real  and  imaginary  parts  of  o’, ,  and,  therefore,  only  one  of  the  roots  was  considered  in  the 
asymptotic  calculation.  The  asymptotic  stress  components  in  Eq.  (24)  can  be  rewritten  for  the 
complex  root  as 

o,.  =f;. /  =  1,A,6 

costo-Im(o,)  sinco}+/  {Re(o.)  sinco  +Im(o.)  cosco}], 

where  co  =  Im(?ij)  ln|r|| .  With  the  complex  number  of  K . ,  Eq.  (30)  becomes 

o .  =  F.  +  Fj  r|'^®*^'^"*{Re(^ )  cosco  -  Im(^ )  sin  co} 

+  Kj  r|’^°^’'‘^“*{Re(o\ )  sin  co  +  Im(o\  )  cosco},  /  =  1,  A  ,  6, 

where  and  are  pure  real  numbers.  Therefore,  the  complex  roots  also  led  to  the  two-term 

asymptotic  solution. 

Figure  38  shows  distributions  of  amplitudes  of  real  and  imaginary  parts  of  the  singular 
stress  components  in  the  global  xyz  -axes  as  functions  of  \\f  for  the  complex  root.  The  three 
interlaminar  stress  components  are  nearly  constant  in  the  real  part,  while  no  constant  component 
is  found  in  the  imaginary  part.  The  traction  continuity  boundary  conditions  were  satisfied  at  the 
ply  interface  (\\f  =  0 )  in  both  parts. 
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(a)  Re(o,.)  (b)  Im(o,.). 


Figure  38.  Distribution  of  Magnitudes  of  Singuiar  Stress  Components  with  a  Compiex 
Root,  =  0.9646  +  0.0968  i ,  for  (|)  =  60° . 


Figures  39  and  40  show  stress  distribution  by  asymptotic  and  numerical  analyses  in  the 
global  xyz  -axes  and  the  local  rst  -axes.  The  comparisons  in  both  figures  were  made  with  the 
same  multiplicative  coefficients,  Ki  =  -2  and  K2  =  0.2  corresponding  to  the  first  and  the  second 
power  of  singularities.  The  additives  in  Figure  39  were  Fi  =  -0.72,  F2  =  0.45,  F3  =  -0.23, 

F4  =  -0.04,  F5  =  0.4,  and  F(,  =  0,  while  in  Figure  40  were  F2  =  0.44,  F3  =  -0.95, 

F4  =  -0.03,  and  Fi  =  F5  =  Fe  =  0. 

In  Figure  39,  with  the  range  of  x/h>0.0\,  excellent  agreement  was  observed  between 
the  full- field  numerical  solution  and  two-term  asymptotic  one  for  all  stresses  with  one  set  of  the 
multiplicative  factors  and  the  additives.  One  exception  was  the  interlaminar  normal  stress 
which  showed  a  rapid  increase  in  the  asymptotic  solution  as  the  point  moved  away  from  the 
singular  point.  All  stress  components  grew  infinitely  on  approaching  the  singular  point  along  the 
interface  line  (  z*  =  ±0 )  in  the  asymptotic  solution,  but  yielded  finite  values  in  the  numerical 
solution. 
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(d)  Interlaminar  shear  stress  (CJ  ) 


(f)  In-plane  shear  stress  (CJ  ^ ) 


Figure  39.  Stress  Distribution  in  Giobai  xjz-Axes  by  Asymptotic  and  Numericai  Anaiyses  near 
the  Singuiar  Point  Aiong  Straight  Lines  Paraiiei  to  the  interface  between  Two  Piies. 
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(c)  Interlaminar  normal  stress  (G^^) 
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(e)  Interlaminar  shear  stress  (G 
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(b)  In -plane  normal  stress  (G  ) 
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(d)  Interlaminar  shear  stress  (G  ) 


(f)  In-plane  shear  stress  (G 


Figure  40.  Stress  Distribution  in  Locai  rst  -Axes  with  ({)  =  30^  by  Asymptotic  and  Numericai 
Anaiyses  Aiong  Contours  of  Constant  Radii  with  the  Center  at  the  Singuiar  Point. 
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In  Figure  40  three  stress  components  (ct  ,  a  „  and  a  „)  from  the  asymptotic  analyses 
were  identically  zero  at  the  free  edge  (\[t  =  ±7t/2+(|) ),  which  satisfied  the  traction- free  boundary 
condition  in  Eq.  (19),  while  the  numerical  solution  yielded  nearly  zero  stresses.  Excellent 
agreement  was  also  observed  between  the  full- field  numerical  solution  and  two -term  asymptotic 
one  for  all  stresses  along  the  three  contours.  One  exception  was  again  the  interlaminar  normal 
stress  ( ct  „ ),  which  showed  more  discrepancy  as  the  radius  of  the  contour  increased. 

4)  =45° 

Two  roots  found  by  the  asymptotic  analysis  were  =  0.8543  and  =  1.1770 .  Both 
roots  were  significantly  influential  to  the  asymptotic  behavior  because  of  the  weak  singularity. 
Figure  41  shows  distributions  of  amplitudes  of  the  singular  stress  components  in  the  global  xyz  - 
axes  as  functions  of  \|/  for  the  two  roots.  None  of  the  stress  components  is  constant  for  both 
roots.  The  traction  continuity  boundary  conditions  were  satisfied  at  the  ply  interface  (\i/  =  0 ). 

Figures  42  and  43  show  stress  distribution  by  asymptotic  and  numerical  analyses  in  the 
global  xyz  -  axes  and  the  local  rst  -axes.  The  comparisons  in  both  figures  were  made  with  the 


Figure  41.  Distribution  of  Magnitudes  of  Singular  Stress  Components  with  Two  Roots  for  (|)  =  45° . 
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(d)  Interlaminar  shear  stress  (CJ  ) 


(e)  Interlaminar  shear  stress  (CJ  ) 


(f)  In-plane  shear  stress  (CJ  ^ ) 


Figure  42.  Stress  Distribution  in  Giobai  xjz-Axes  by  Asymptotic  and  Numericai  Anaiyses  near 
the  Singuiar  Point  Aiong  Straight  Lines  Paraiiei  to  the  interface  between  Two  Piies. 
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(f)  In-plane  shear  stress  (C 


Figure  43.  Stress  Distribution  in  Locai  rst  -Axes  with  \|/  =  45^  by  Asymptotic  and  Numericai 
Anaiyses  Aiong  Contours  of  Constant  Radii  with  the  Center  at  the  Singuiar  Point. 


72 


same  multiplicative  coefficients,  =  0.45  and  =  0.9  ,  corresponding  to  the  first  and  the 
second  power  of  singularities.  Two  multiplicative  coefficients  were  in  the  same  order  of 
magnitude,  which  indicated  that  both  roots  contributed  similarly  to  the  comparison  of  the  two 
solutions.  The  additives  in  Figure  42  were  F\  =  -1.25,  F2  =  -0.1,  F3  =  -1.25,  =  1.25  and 

F4  =  Fg  =  0 ,  and  in  Figure  43  were  F2  =  -0.1,  F3  =  -2.5  and  F^  —  F^-  F^  —F^—0. 

In  Figure  42,  with  the  range  of  x//i  >  0.01 ,  excellent  agreement  was  observed  between 
the  full- field  numerical  solution  and  two-term  asymptotic  one  for  all  stresses  with  one  set  of  the 
multiplicative  factors  and  the  additives.  One  exception  was  the  interlaminar  normal  stress  (<7^.^), 
which  showed  a  rapid  increase  in  the  asymptotic  solution  as  the  point  moved  away  from  the 
singular  point.  All  stress  components  grew  infinitely  on  approaching  the  singular  point  along  the 
interface  line  ( z*  =  ±0 )  in  the  asymptotic  solution,  but  yielded  finite  values  in  the  numerical 
solution. 

In  Figure  43  three  stress  components  and  from  the  asymptotic  analyses 

were  identically  zero  at  the  free  edge  (\|/  =  ±7t/2-l-(|) ),  which  satisfied  the  traction- free  boundary 
condition  in  Eq.  (19),  while  the  numerical  solution  yielded  nearly  zero  stresses.  Excellert 
agreement  was  also  observed  between  the  full- field  numerical  solution  and  two-term  asymptotic 
one  for  all  stresses  along  the  three  contours.  One  exception  was  again  the  interlaminar  normal 
stress  (o,,),  which  showed  more  discrepancy  as  the  radius  of  the  contour  increased. 

^  =60° 

Two  roots  found  by  the  asymptotic  analysis  were  =  0.7562  and  X2  =  1.4538  .  Both 
roots  were  significantly  influential  to  the  asymptotic  behavior  because  of  the  weak  singularity. 
Figure  44  shows  distributions  of  amplitudes  of  the  singular  stress  components  in  the  global 
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Figure  44.  Distribution  of  Magnitudes  of  Singular  Stress  Components  with  Two  Roots  for  (|)  =  60° . 


xyz  -axes  as  functions  of  \\f  for  the  two  roots.  None  of  the  stress  components  is  constant  for 
both  roots.  The  traction  continuity  boundary  conditions  were  satisfied  at  the  ply  interface 
(¥=0). 

Figures  45  and  46  show  stress  distribution  by  asymptotic  and  numerical  analyses  in  the 
global  xyz  -axes  and  the  local  rst  -axes.  The  comparisons  in  both  figures  were  made  with  the 

same  multiplicative  coefficients,  =  0.1 1  and  A' 2  =0.1,  corresponding  to  the  first  and  the 
second  power  of  singularities.  Two  multiplicative  coefficients  were  in  the  same  order  of 
magnitude,  which  indicated  that  both  roots  contributed  similarly  to  the  comparison  of  the  two 
solutions.  The  additives  in  Figure  45  were  Fi  =  -0.72,  F2  =  0.45  ,  F3  =  -0.23,  F4  =  -0.04, 

Fj  =  0.4  and  F^  =  0 ,  and  in  Figure  46  were  F^  =  0.44 ,  F3  =  -0.95,  F4  =  -0.03  and 
F,=F,=F,  =0. 

In  Figure  45,  with  the  range  of  x/h  >  O.Ol,  excellent  agreement  was  observed  between 
the  full- field  numerical  solution  and  two-term  asymptotic  one  for  all  stresses  with  one  set  of  the 
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Figure  45.  Stress  Distribution  in  Giobai  xjz-Axes  by  Asymptotic  and  Numericai  Anaiyses  near 
the  Singuiar  Point  Aiong  Straight  Lines  Paraiiei  to  the  interface  between  Two  Piies. 
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Figure  46.  Stress  Distribution  in  Locai  rst  -Axes  with  ({)  =  60^  by  Asymptotic  and  Numericai 
Anaiyses  Aiong  Contours  of  Constant  Radii  with  the  Center  at  the  Singuiar  Point. 
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multiplicative  factors  and  the  additives.  One  exception  was  the  interlaminar  normal  stress 
which  showed  a  rapid  increase  in  the  asymptotic  solutbn  as  the  point  moved  away  from  the 
singular  point.  All  stress  eomponents  grew  infinitely  on  approaching  the  singular  point  along  the 
interface  line  (  z*  =  ±0 )  in  the  asymptotic  solution  but  yielded  finite  values  in  the  numerical 
solution. 

In  Figure  46  three  stress  components  ( and  )  from  the  asymptotie  analyses 
were  identically  zero  at  the  free  edge  (\|/  =  ±7t/2+(|) ),  which  satisfied  the  traction- free  boundary 
eondition  in  Eq.  (19),  while  the  numerical  solution  yielded  nearly  zero  stresses.  Exeellent 
agreement  was  also  observed  between  the  full- field  numerical  solution  and  two-term  asymptotic 
one  for  all  stresses  along  the  three  contours.  One  exception  was  again  the  interlaminar  normal 
stress  (o,,),  which  showed  more  diserepaney  as  the  radius  of  the  contour  increased. 

6.4  SUMMARY  AND  RECOMMENDATIONS 

Asymptotic  and  numerical  analyses  were  performed  to  provide  an  accurate  stress  field  of 
laminated  composites  with  slanted  free  edges  in  the  vicinity  of  the  ply  interface  and  the  free 
edge.  The  analyses  were  performed  with  various  slanted  angles.  The  numerical  analysis  was 
developed  based  on  the  B-spline  approximation  method  of  displaeement  funetions.  Stress 
components  by  the  numerical  method  were  obtained  for  a  laminated  composite  with  the  lay-up  of 
[45/-45]s  AS4/3501-6  laminate  under  a  uniaxial  loading. 

The  behavior  of  power  of  singularity  was  highly  dependent  on  the  slanted  angles  of  the 
free  edge.  The  power  of  singularities  was  found  with  two  real  numbers  at  0°  <  (|)  <  2°  and 
38°  <  (|)  <  90° ,  and  with  complex  conjugated  numbers  at  3°  <  (|)  <  37° .  At  every  slanted  angle, 
excellent  agreement  was  observed  for  the  stress  distributions  to  a  distance  of  approximately 
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one-half- ply  thickness  from  the  singular  point  between  the  full- field  numerical  solution  and  the 
two-term  singular  asymptotic  solution  with  appropriate  multiplicative  coefficients  and  constant 
additives. 
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7.  MECHANICAL  MODELING  AND  PREDICTION  OE  BULK  PROPERTIES  OE 

OPEN-CELL  CARBON  FOAM 

Carbon  foam  is  emerging  as  an  ultralightweight  and  efficient  thermal  management 
material  for  many  multifunctional  applications,  such  as  high-density  electronics,  hybrid  diesel- 
electric  vehicles,  communication  satellites,  and  advanced  aircraft  [37],  Carbon  foam  was  shown 
to  demonstrate  numerous  unique  properties  that  make  it  an  attractive  material  for  the  use  of  the 
low-cost,  lightweight,  insulating,  energy- absorbing  structural  components  [38], 

The  carbon  foam  is  expected  to  possess  isotropic  properties  macroscopically.  However,  a 
microstructure  of  an  open- cell  foam,  based  on  the  minimum  surface  energy  during  the  foaming 
process,  possesses  a  tetrahedral  structure  of  the  foam  ligaments  oriented  approximately  109 
degrees  with  each  other  [39,40],  as  in  Figure  47.  Preliminary  investigation  reveals  that  the 
graphitic  alignment  of  the  cell  ligaments  varies  along  its  longitudinal  (axial)  direction. 

Processing  parameters,  such  as  temperature,  pressure,  etc.,  control  the  porosity  and  the  graphitic 
alignment  of  the  carbon  foam,  which  in  turn  determine  its  geometries  and  material  properties. 
With  an  ongoing  research  effort  to  investigate  an  appropriate  microstructural  characterization 
technique  to  correlate  the  foam  microstructure  with  the  processing  parameters,  the  micro- 
structural  geometry  and  material  properties  of  the  foam,  including  mechanical  elastic  moduli, 
Poisson’s  ratio,  thermal  conductivity,  etc.,  will  be  used  for  the  mechanical  and  thermal  analysis. 

Because  of  the  tetrahedral  cell  microstructure  of  the  carbon  foam,  the  macroscopic 
properties,  such  as  foam  moduli  and  strengths,  are  critically  influenced  by  the  deformation 
characteristics  of  the  cell  ligaments.  Further,  as  stated  earlier,  the  cross-sectional  area  and  the 
graphitic  alignment  (hence  the  directional  material  properties)  of  the  cell  ligaments  vary  along 
the  ligament  length.  The  tetrahedral  cell  microstructure  consists  of  four  ligaments  as  a  frame 
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Figure  47.  Microstructure  of  Carbon  Foam. 

structure,  and  the  four  ligaments  connected  to  a  common  node  are  oriented  approximately  109 
degrees  to  one  another  in  3-D  space.  The  cross-sectional  area  and  the  material  properties  vary 
along  the  longitudinal  directions  of  the  ligaments.  Because  of  the  complex  geometry  and 
anisotropic  material  properties,  it  is  appropriate  to  perform  the  analysis  numerically  to  obtain 
accurate  displacement  and  stress  field  solutions. 

The  numerical  model  should  be  able  to  analyze  the  deformation  of  the  ligaments  that  are 
connected  and  oriented  in  the  3-D  space  under  arbitrary  loading  conditions.  The  analysis  should 
be  able  to  predict  longitudinal  and  transverse  displacements  as  well  as  rotations,  and  to  calculate 
the  stress  and  strain  distributions  along  the  ligaments.  To  achieve  the  above  goals,  the  geometry 
of  cell  ligaments  in  the  carbon  foam  needs  to  be  characterized  based  on  the  bubble -forming 
process.  With  the  geometry  known,  the  stress  analysis  can  be  performed  with  the  microstructure 
of  the  carbon  foam  under  a  certain  boundary  condition.  By  applying  appropriate  boundary 
conditions,  we  can  then  calculate  effective  elastic  Young’s  moduli  and  Poisson’s  ratio  of  the 
foam  in  the  macroscopic  level. 
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7.1  MODELING  OF  CARBON  FOAM 


7.1.1  Characterization  of  Geometry  of  Unit  Cell  of  Carbon  Foam 

Because  of  the  randomness  and  complexity  of  the  microstructure  of  the  carbon 
foam,  it  is  difficult  to  consider  whole  individual  cells  in  the  foam.  Instead,  an  assumption  is 
made  that  the  behavior  of  the  material  is  represented  with  a  unit  cell  of  the  foam  containing  four 
ligaments.  The  microstructural  characterization  of  the  representative  unit- cell  ligaments  will 
then  be  correlated  with  the  macroscopic  bulk  properties  in  a  statistical  means.  To  obtain  a 
representative  response  of  the  unit  cell,  the  geometry  and  pertinent  material  properties  of  the  unit 
cell  are  based  on  the  collection  of  data  on  the  size,  shape,  and  property  variation  of  the  foam 
ligaments. 

To  analyze  the  materials  behavior  of  the  foam,  we  need  to  generate  the  geometry 
of  the  unit  cell  of  the  foam  in  an  appropriate  manner.  The  unit  cell  can  be  generated  by  making 
use  of  symmetry  that  exists  in  the  bubble- forming  process.  Surface  energy  of  nucleation  and 
expansion  of  the  bubbles  during  the  foam- forming  process  are  minimal  when  the  centers  of  the 
bubbles  coincide  with  the  four  vertices  of  a  tetrahedron.  The  vertices  of  the  tetrahedron  are 
defined  by  first  considering  a  cube.  For  the  sake  of  simplicity  of  defining  the  unit  cell,  the 
dimension  of  the  cube  is  taken  as  2ax2ax2a  in  the  x-,  y-  and  z-directions,  with  its  origin  (point 
0)  located  at  the  center  of  the  cube,  as  in  Figure  48.  The  vertices  of  the  tetrahedron  confined  in 
the  cube  are  the  four  corner  points  of  the  cube  that  are  located  diagonally  with  each  other  on  the 
faces  of  the  cube  (e.g.,  points  2,  4,  5,  7  in  Figure  48).  Connecting  the  four  comer  points  then 
generates  a  tetrahedron,  whose  volume  is  =So^ /3.  This  tetrahedron  can  be  divided  into 

four  equal  subregions  containing  the  origin  of  the  tetrahedron.  The  vertices  of  the  subregions  are 
shown  in  Figure  48  as  0-2-7-4,  0-7- 5-4,  0-4- 5- 2,  and  0- 2-7-5,  respectively. 
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Figure  48.  Points  and  Lines  in  a  Cube  and  a  Tetrahedron  for  Generating  a  Unit  Celi  of  Carbon 
Foam. 

The  next  step  is  to  generate  four  spheres  growing  from  the  four  vertices  (the  four 
corner  points)  of  the  tetrahedron.  The  spheres  represent  bubbles  that  are  produced  during  the 
foaming  process.  The  radii  of  the  spheres  will  determine  the  porosity  of  the  unit  cell  of  the  foam. 
Subtracting  the  volume  of  the  bubbles  (spheres)  from  that  of  the  tetrahedron  creates  the  unit  cell 
of  the  carbon  foam,  as  Figure  49  shows.  The  porosity  ( (|) )  of  the  foam  is  thus  calculated  by 

(j)  -i—V^^itlVtetra  5  where  is  the  volume  of  the  unit  cell.  Figure  50  shows  the  unit  cells  with 
three  different  levels  of  porosity. 

For  convenience,  local  coordinate  systems,  whose  x-directions  are  parallel  to  the 
longitudinal  directions  (axes)  of  the  ligaments,  are  defined  by  using  four  lines  that  connect  the 
points  (0-15,  0-16,  0-17,  and  0-18)  in  Figure  48.  The  local  coordinate  systems  are  useful  because 
the  ligaments  possess  material  symmetry  along  its  longitudinal  axes  due  to  the  graphitic 
alignment. 

For  the  numerical  analysis  using  an  FEM,  the  ligaments  of  the  foam  generated 
above  are  discretized  into  FEs  along  the  local  longitudinal  directions. 
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Figure  49.  A  Tetrahedron  and  Spheres  to  Generate  a  Unit  Cell  of  Carbon  Foam. 


(a)  (|)  =  10%  (b)  (|)  =  78%  (c)  (|)  =  98% 

Figure  50.  Unit  Celis  of  Carbon  Foam  with  Various  Porosities. 


The  coordinates  of  the  cross-sectional  area  as  well  as  the  nodal  points  along  the  longitudinal 
directions  are  calculated  by  the  following  method:  a  shaded  volume  in  Figure  48  is  selected  as 
one  of  the  four  subtetrahedra  and  plotted  in  Figure  51.  The  selected  subtetrahedron  can  be 
subdivided  further  into  three  sub -subtetrahedra  (0-15-2-7,  0-15-7-4,  and  0-15-4-2)  with  respect 
to  symmetry  planes. 

As  the  bubble  grows  from  the  three  vertices  of  the  subtetrahedron,  increasing  the 
porosity,  the  unit  cell  undergoes  three  different  shapes.  Figure  52  shows  the  three  different 
cross-sectional  slices  in  plane  2-4-7  during  the  bubble -forming  process.  When  the  porosity  of 
the  carbon  foam  is  78  percent,  the  bubbles  contact  with  each  other.  For  the  open-cell  carbon 
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Figure  51.  One  of  Four  Subtetrahedra  Considering  Symmetry  in  Figure  48. 


Figure  52.  Cross-Sectional  Slices  During  Bubble-Forming  Process  (Ris  the  Radius  of  Spheres). 


foams,  whose  typical  porosities  are  higher  than  80  percent,  the  bubbles  coalesce  with  each  other, 
as  shown  in  Figure  52(c). 

Figure  53  shows  the  intersection  of  the  subtetrahedron  and  the  bubble  spheres. 
The  subtetrahedron  subtracted  from  the  bubbles  defines  the  ligament  of  the  carbon  foam.  As 
Figure  53  shows,  points  A,  B  and  a,  b  are  the  intersections  of  the  half  of  the  sub -subtetrahedron 
(0-15-9-7)  and  a  bubble  sphere  with  the  radius  of  R.  The  coordinates  of  these  points  are  as 
follows: 
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The  subtracted  sub -subtetrahedron  in  Figure  53  is  further  discretized  into  the  finite  volumes 


along  the  ligament  axis  (  P^P2  ).  Discretized  points  (Qi,Q^  ,A  )  on  the  ligament  axis  form 


discretization  planes  perpendicular  to  the  ligament  axis,  which  intersect  with  the  lines  (P^P^  and 


P^Pf^ )  on  the  ligament  surface  at  points  5j  ,S^  ,A  and  Vj  ,F^,A  ,  respectively,  as 

shown  in  Figure  54.  The  coordinates  of  the  points  ( 2, ,  and  F, )  are  obtained  by  using 
geometric  arithmetic  as  follows: 
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Figure  54.  A  Sub-subtetrahedron  Discretized  into  Finite  Volumes  Along  a  Ligament  Axis. 


The  coordinates  of  the  points  in  other  sub-subtetrahedra  are  obtained  by  reflection 
with  respect  to  symmetry  planes.  For  example,  a  point  5-  on  a  plane  (0-15-12)  is  the  reflection 


of  the  point  5,  on  a  plane  (0-15-9)  with  respect  to  a  symmetry  plane  (0-15-7),  as  can  be  inferred 
from  Figure  51.  The  formula  for  the  reflection  in  a  plane  with  equation  ax+by  +  cz  +  d  =  0  is 


(x,y)  a  ^ - - - j[M{x„,y„,Zo)-{2ad2bd,2cd)], 

a  +b  +c 


where  M  is  a  matrix  given  by 
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and  the  formula  for  the  plane  going  through  (x^,y^,Zg),  ( 


Xj,y,,Zi)and(x2,y2,Z2)i^  given  by 
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Once  the  coordinates  of  the  points  in  Figure  54  (c)  are  obtained  by  using  Eqs.  (33) 
through  (36),  the  ligament  geometry  of  the  unit  cell  as  well  as  the  material  property  distribution 
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along  the  ligaments  can  be  generated  with  the  nodal  points  ( 2, )  with  the  cross-sectional 
information  (  V' ,  S" , V" ). 

The  above  procedure  of  defining  the  model  geometry  can  be  used  in  either 
developing  a  simplistic  FE  model,  assuming  the  ligaments  behave  like  beam  elements,  or 
generating  3-D  FE  meshes  with  commercial  FE  model  packages,  such  as  ANSYS.  Figure  55 
shows  3-D  FE  meshes  generated  by  ANSYS  and  equivalent  1-D  beam  meshes.  The  solution  of 
the  beam  FEA  requires  further  development  of  the  beam  model  and  thus  is  not  included  in  this 
section.  However,  the  stress  analysis  of  the  unit  cell  using  the  3-D  meshes  is  presented  in  this 
section. 


(a)  3-D  mesh  (b)  1-D  mesh 

Figure  55.  3-D  Mesh  versus  1-D  Beam  Mesh  of  the  Unit  Ceil  of  Carbon  Foam. 


7.1.2  Determination  of  Boundary  Condition 

The  FEA  requires  appropriate  loading  and  boundary  conditions  to  be  applied  on 
the  unit  cell  of  the  foam.  Eor  the  analysis  of  the  unit  cell,  we  need  to  correlate  the  overall 
loading/boundary  conditions  (OBC)  with  the  ligament  boundary  conditions  (LBC).  The  LBC 
varies  with  the  location,  size  and  orientation  of  the  unit  cell  in  the  carbon  foam  under  an  OBC. 
However,  to  correlate  LBC  with  OBC,  the  macroscopic  (bulk)  material  symmetry  of  the  foam 
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needs  to  be  known.  The  measured  macroscopic  (bulk)  properties  of  several  carbon  foams 
indicate  that  the  bulk  properties  of  carbon  foam  are,  in  general,  isotropic  [41,42],  The  bulk 
material  isotropy  is  thus  assumed  in  this  study. 

The  LBC  can  be  determined  by  generating  an  imaginary  tetrahedron  inside  which 
a  unit  cell  of  the  foam  is  to  be  located,  as  in  Figure  56.  The  OBCs  are  then  applied  to  the 
surfaces  of  the  tetrahedron.  Under  the  OBC,  the  bulk  tetrahedron  is  assumed  to  be  isotropic,  and 
the  displacements  in  the  tetrahedron  are  calculated  point  by  point.  The  displacements  are 
collected  for  the  points  that  coincide  with  the  tips  of  the  ligament.  The  FEA  is  then  run  by 
assigning  the  corresponding  displacement  boundary  conditions  to  the  unit  cell  that  was  generated 
in  the  previous  section. 

The  FEA  is  repeated  with  the  variation  of  size  and  orientation,  as  the  unit  cells  are 
distributed  randomly  in  the  bulk  of  carbon  foam.  As  the  selection  of  the  input  parameters  (size, 
orientation)  increases,  the  computational  time  for  numerical  analysis  increases  rapidly,  and 
therefore,  we  need  the  1-D  beam  analysis  instead  of  the  3-D  analysis.  A  statistical  approach 
should  be  incorporated  for  the  random  cell  size  and  its  orientations. 

7.2  NUMERICAL  CALCULATION 

7.2.1  Prediction  Method  of  Effective  Bulk  Moduli  of  Foam 

To  assess  the  validity  of  the  model,  effective  bulk  moduli  were  predicted  for 
polyurethane  foam  as  well  as  the  carbon  foam,  based  on  the  ligament  material  properties. 
Although  the  material  orthotropy  of  the  ligaments  can  be  included  in  the  effective  properties 
prediction,  at  present  for  convenience,  ligaments  are  assumed  as  isotropic  material.  Densities  of 
the  ligaments  (  in  the  polyurethane,  carbonized  and  graphitized  foams  were  taken  as  1.2, 


Figure  56.  An  Imaginary  Tetrahedron  inside  which  a  Unit  Cell  of  the  Foam  is  Located. 

1.8,  and  2.2  g/cm^,  respectively.  The  foam  bulk  density  ( can  then  be  related  with  the 
material  density  of  the  ligament  as 

P bulk  ~  P ligament  (1  •  (37) 

Since  the  size  of  the  bubbles  determined  the  porosity,  the  foam  density  can  be  related  with  the 
bubble  size  (radius)  as  well,  as  plotted  in  Figure  57. 

The  FEA  was  performed  with  ANSYS.  The  FE  meshes  of  the  unit  cell  were 
generated  with  unit  length  of  cube,  tetrahedron  and  bubble  spheres  with  radius  determined  by 
Figure  57.  The  unit  cell  was  loaded  under  a  constant  uniaxial  strain  in  a  loading  direction,  and 
constrained  to  prevent  deformation  in  the  other  two  perpendicular  directions.  The  uniform 
uniaxial  strain  boundary  condition  was  achieved  by  applying  a  prescribed  displacement 
boundary  condition  on  nodal  points  on  the  tip  surfaces  of  the  ligaments,  where  the  tip  surfaces 
coincided  with  the  surfaces  of  the  tetrahedron,  as  explained  earlier.  The  values  of  the  prescribed 
displacements  in  the  loading  direction  varied  with  the  location  of  the  nodal  points,  whereas  those 
in  other  directions  perpendicular  to  the  loading  direction  were  set  to  zero.  For  example,  if  the 
loading  was  applied  in  the  x-direction,  the  displacement  values  were 
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Figure  57.  Relation  between  Bulk  Density,  Porosity  and  Bubble  Size  for  Open-Cell  Foams. 


d  =  J  =  0  ,  d 

X  /  X  ’  y  ’  z 


=  0  , 


(38) 


where  x  and  e "  are  the  nodal  coordinate  and  the  applied  constant  strain,  respectively. 

Effective  stress  components  ( ^ )  can  be  calculated  by  taking  volumetric  averages 
of  resultant  stress  field  over  the  volume  of  the  tetrahedron,  such  as 


a,V 


CT,  =■ 


a,V 


E  . 

unit 


(39) 


where  g/  and  v"  are  the  elemental  stress  and  volume  of  the  unit  cell,  respectively.  According  to 

the  3-D  constitutive  law  of  isotropic  material,  the  effective  moduli  ( E )  and  Poisson’s  ratio  (v”) 
of  the  unit  cell  were  calculated  as  follows: 
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^[g,-v(g.+g,)] 


V  =■ 


Ct.+G, 


(40) 
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where  a ,  is  the  volumetric  average  of  stress  components  in  the  loading  direction,  and  a  j  and 

are  those  in  the  other  two  directions  perpendicular  to  the  loading  direction. 

The  carbon  foam  is  filled  with  the  unit  cells  that  are  randomly  oriented  in  3-D 
space  and  connected  with  each  other  at  the  tips  of  the  ligaments.  The  random  nature  of 
orientation  leads  to  the  isotropic  properties.  Therefore,  collection  of  the  randomly  oriented  unit 
cells  needed  to  be  considered  in  the  calculation  of  the  effective  bulk  moduli  of  the  carbon  foam. 
Instead  of  the  rigorous  random  generation  of  the  orientation,  however,  different  orientation 
schemes  were  employed  in  this  study  in  a  less  rigorous  manner  using  mutually  orthogonal 
coordinate  systems:  Cartesian  and  spherical  coordinate  systems,  as  shown  in  Figure  58. 

When  the  Cartesian  coordinate  system  was  used  to  rotate  the  unit  cell  with  respect 
to  X-,  y-  and  z-axes,  as  shown  in  Figure  58(a),  only  two  axes  among  the  three  axes  were 
considered,  since  the  orientation  with  respect  to  the  axis  parallel  to  the  uniaxial  loading  did  not 
alter  the  effective  moduli  of  the  unit  cell.  For  example,  when  the  load  was  applied  to  the 
direction  parallel  to  the  x-axis,  only  two  rotations  with  respect  to  y-  and  z-axes  were  considered 
in  the  moduli  calculation.  Meanwhile,  with  the  spherical  coordinate  system,  two  angles,  zenith 
(0  )  and  azimuth  ((|) )  angles,  were  used  to  rotate  the  unit  cell,  as  shown  in  Figure  58(b).  Three 
different  rotating  schemes  tried  in  these  coordinate  systems  were 

1.  Scheme  I:  first  rotation  in  y-axis  ( 0  ^  )and  the  second  rotation  in  z-axis  ( 0^ ). 

2.  Scheme  II:  first  rotation  in  z-axis  ( 0^ )  and  the  second  rotation  in  y-axis  (0  ). 

3.  Scheme  III:  first  rotation  in  zenith  (0  )  and  the  second  rotation  in  azimuth  ((|) ). 

Young’s  moduli  of  the  unit  cells  were  calculated  by  using  Eq.  (40)  for  an 
isotropic  foam  with  the  changes  of  angles  in  Schemes  I,  II  and  III,  and  the  results  were  plotted  in 
Figures  59,  60,  and  61,  respectively.  The  bulk  Young’s  modulus  and  Poisson’s  ratio  of  the 
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(a)  Cartesian  coordinate 


(b)  Spherical  coordinate 


Figure  58.  Coordinate  Systems  for  Rotating  Unit  Celi  for  Simulating  Random  Orientation. 


(a)  Surface  plot 


(b)  Contour  plot 


Figure  59.  Variation  of  Young’s  Moduii  of  Unit  Celis  of  Isotropic  Foam  on  the  Changes  of  Angles 
in  Cartesian  Coordinate  Axes,  with  the  First  Rotation  in  the  y-Axis  and  the  Second 
Rotation  in  the  z-Axis. 


(a)  Surface  plot  (b)  Contour  plot 


Figure  60.  Variation  of  Young’s  Moduii  of  Unit  Celis  of  Isotropic  Foam  on  the  Changes  of  Angles 
in  Cartesian  Coordinate  Axes,  with  the  First  Rotation  in  the  z-Axis  and  the  Second 
Rotation  in  the  3/^Axis. 
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(a)  Surface  plot  (b)  Contour  plot 

Figure  61.  Variation  of  Young’s  Moduii  of  Unit  Celis  of  Isotropic  Foam  on  the  Change  of  Angies  in 
Sphericai  Coordinate  Axes. 

material  were  4.5  GPa  and  0.33,  respectively.  Both  mesh  and  contour  plots  were  drawn  for  the 
orientation  angles  from  0°  to  360°.  Average  Young’s  moduli  and  Poisson’s  ratios  are  listed  in 
Table  8.  The  figures  show  periodic  patterns  of  the  moduli  in  both  rotating  angles  either  in  a 
symmetric  way  or  an  antisymmetric  way.  The  periodic  variations  observed  in  each  scheme  were 
as  follows: 

1.  Scheme  I: 

£(0,,0J  =  £(e,,0  +18O°)  =  £(0,  +18O°,18O°-0J  =  £(0  +18O°,36O°-0  A  (41) 

2.  Scheme  II: 

£(0  „  ,0 , )  =  £(0  + 1 80°, 0  J  =  £(1 80°  -0  „  ,  0  + 1 80°)  =  £(360°  -0  „  ,  0  + 1 80°)  (42) 

y  z  y  z  y  z  y  z 

3.  Scheme  III: 

£(0  ,^)  =  £(0  + 180°,  4))  =  £(1 80°  -0 , 4)  + 1 80°)  =  £(360°  -0 , 4)  + 1 80°)  (43) 

Because  of  the  periodicity,  to  completely  fill  the  isotropic  bulk  foam  with  the  unit 
cells,  it  was  necessary  to  rotate  the  unit  cell  only  from  0°  to  180°  in  both  axes.  Since  the  moduli 
were  calculated  by  taking  the  average  over  all  orientations,  it  was  not  sensitive  on  the  selection 
of  the  coordinate  system.  In  this  study  the  spherical  coordinate  system  was  chosen  as  the 
rotating  method  in  the  moduli  calculation. 
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Table  8 

Predicted  Effective  Young’s  Moduii  and  Poisson’s  Ratio  of  an  isotropic  Foam 


E  (MPa) 

V 

Scheme  1 

251.9 

0.37 

Scheme  II 

241.9 

0.39 

Scheme  III 

251.9 

0.37 

7.2.2  Results:  Polyurethane  Foam 

Effective  bulk  properties  of  the  polyurethane  foam  were  calculated  with  the 
present  method  at  various  densities.  The  Young’s  modulus  and  Poisson’s  ratio  of  the 
polyurethane  as  the  ligament  properties  were  1.8  GPa  and  0.33,  respectively.  The  property 
predictions  by  the  present  method  were  compared  with  those  both  in  a  handbook  [37]  and  an 
existing  formula  suggested  by  Warren  and  Kraynik  [39],  as  shown  in  Figure  62.  The  latter  was 
formulated  based  on  a  similar  unit-cell  model  as  the  present  one,  but  limited  to  the  ligaments  that 
possess  the  isotropic  material  properties  and  uniform  cross  sections,  such  as  circular,  triangular, 
and  plateau  border  cross  sections.  The  plateau  border  cross  section  corresponds  to  the  space 
between  three  identical,  mutually  tangent  circles.  Meanwhile,  the  handbook  provides  only  shear 
moduli  at  various  ranges  of  densities.  Therefore,  the  Young’s  moduli  were  back  calculated  from 
the  shear  moduli,  along  with  an  assumed  Poisson’s  ratio  of  0.33.  Note  that  since  the  radius  of 

the  bubbles  for  the  open-cell  foams  was  greater  thanV2a,  the  bulk  density  of  the  open-cell 
foams  should  be  less  than  0.266  g/cm^. 

For  low-density  foams,  the  bending  compliance  of  the  ligaments  is  much  greater 
than  the  axial  compliance.  Therefore  the  bulk  Young’s  modulus  is  dominated  by  the  bending 
behavior  of  foam  ligaments.  The  dominance  on  the  bending  behavior  results  in  quadratic 
variation  of  the  modulus  with  the  density  increment,  as  was  also  concluded  in  the  previous 
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Figure  62.  Prediction  of  Buik  Young’s  Moduii  of  Polyurethane  Foam  at  Various  Densities. 


prediction  method  by  Warren  and  Kraynik.  As  lying  on  upper  bound,  the  present  results  agreed 
well  with  both  the  existing  formula  and  the  measured  data. 

Figure  63  shows  predicted  Poisson’s  ratios  of  the  polyurethane  foam  at  various 
densities.  At  the  nearly  zero  foam  density,  both  predictions  of  the  Poisson’s  ratio  approached 
0.5.  While  Warren’s  predictions  decreased  monotonically  with  the  increase  of  the  density,  the 
present  prediction  kept  constant,  or  slightly  increased,  except  at  the  low- density  range  where  the 
Poisson’s  ratio  changed  rapidly.  However,  note  that  the  moduli  by  the  Warren’s  method  were 
calculated  with  the  formula  for  low- density  foams  and  were  extrapolated  to  the  high- density 
range.  Therefore,  we  can  conclude  that  the  Poisson’s  ratio  of  the  polyurethane  foam  was 
dominated  by  the  bending  behavior  at  low-density  range,  but  became  influenced  by  the  axial 
compliance  as  the  density  increased. 

7.2.3  Results:  Carbon  Foam  with  Isotropic  Properties 

Effective  bulk  properties  of  the  carbon  foam  were  calculated  at  various  densities. 
Since  the  material  properties  were  not  well  known  at  the  ligament  level,  the  prediction  was  made 
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Figure  63.  Prediction  of  Buik  Poisson’s  Ratio  of  Poiyurethane  Foam  at  Various  Densities. 

both  on  the  isotropic  and  anisotropic  assumptions  on  the  foam  ligament  properties.  This  section 
first  demonstrated  the  results  based  on  the  isotropic  assumption.  Poisson’s  ratio  was  assumed  to 
be  0.33.  The  density  of  the  carbon  material  was  1.8  g/cm^,  assuming  the  foam  was  carbonized 
but  not  graphitized.  Note  that  since  the  radius  of  the  bubbles  for  the  open-cell  foams  was  greater 

than  42a ,  the  bulk  density  of  the  open-cell  carbon  foams  should  be  less  than  0.398  g/cm?. 

Effective  tensile  modulus  predicted  by  the  present  method  as  well  as  those  by 
Warren  for  the  isotropic  foams  are  shown  in  Figure  64.  The  Young’s  modulus  of  the  foam 
ligament  was  assumed  15.61  GPa,  which  was  the  transverse  Young’s  modulus  of  a  pitch-based 
carbon  fiber,  AS4  [43] .  Both  predictions  resulted  in  quadratic  variation  of  the  moduli  with  the 
density  increment.  The  present  method  for  the  carbon  foam  agreed  well  with  Warren’s 
prediction  with  plateau  border  cross  section  at  a  low-density  region  and  with  the  circular  cross 
section  at  a  high-density  region.  Note  that  for  the  polyurethane  foam,  the  present  method  yielded 
the  upper  bound  of  Warren’s  prediction  at  all  ranges  of  density. 
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Figure  64.  Prediction  of  Buik  Tensile  Moduli  of  Carbon  Foam  at  Various  Densities. 

To  understand  the  influence  of  the  ligament  properties  on  the  bulk  foam 
properties,  the  Young’ s  modulus  of  the  foam  ligament  varied  from  7,  20,  and  25  GPa,  and  the 
corresponding  bulk  moduli  at  various  densities  were  plotted  in  Figure  65.  The  prediction  was 
compared  with  limited  experimental  data  measured  with  MER  carbon  foam  [44].  Figure  65 
shows  that  the  present  method  yielded  a  fairly  reasonable  prediction  of  the  effective  bulk 
properties  of  the  carbon  foam.  Note  that  the  present  method  can  also  be  utilized  to  estimate  the 
ligament  properties  of  the  foam  when  the  bulk  properties  were  measured  by  experiments. 

Figure  66  shows  the  predicted  Poisson’s  ratio  of  the  carbon  foam  at  various 
densities.  As  in  the  case  of  the  polyurethane  foam,  the  present  method  yielded  nearly  constant 
values  unless  the  carbon  foam  was  of  low  density. 

7.2.4  Results:  Carbon  Foam  with  Anisotropic  Properties 

Microstructural  observation  on  the  carbon  foam  revealed  that  the  foam  ligaments 
had  anisotropic  material  properties,  and  the  degree  of  anisotropy  varied  with  foam  processes, 
such  as  carbonization  or  graphitization  The  isotropic  assumption  made  earlier  in  this  study  was 
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Figure  65.  Prediction  of  Buik  Tensile  Moduli  of  Carbon  Foam  at  Various  Densities  with  Several 
Assumed  Ligament  Properties  of  Carbon  Material.  The  present  model  was  compared 
with  experimental  data  [44]. 


Figure  66.  Prediction  of  Bulk  Poisson’s  Ratio  of  Carbon  Foam  at  Various  Densities. 
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dropped  by  varying  the  properties  in  the  longitudinal  and  transverse  directions  parallel  and 
perpendicular  to  the  foam  ligament  axes,  respectively.  With  the  absence  of  well-known  bulk 
properties  of  the  carbon,  several  assumptions  were  made  as  follows:  (1)  the  foam  ligaments 
possessed  transversely  anisotropic  properties;  (2)  in-plane  and  out-of-plane  major  Poisson’s 
ratios  were  assumed  as  v  j2  =Vj3  =  0.33  and  v  23  =  0.5 ,  respectively;  and  (3)  in-plane  ( Gj2  and 

Gjj)  and  out-of- plane  ( G23 )  shear  moduli  were  calculated  by 


<^12  -  - 


2(1+V2i) 


and  G23  = 


2(1+V23) 


(44) 


where  v  21  =v  j2  E^jE^  . 

Effective  bulk  properties  of  the  carbon  foam  predicted  with  the  present  method  at 
various  densities  were  plotted  in  Figures  67  through  71.  Results  with  the  anisotropic  material 
properties  were  compared  with  those  with  isotropic  ones.  The  bulk  moduli  of  the  carbon  foam 
were  predicted:  (1)  when  the  longitudinal  modulus  of  the  ligament  was  fixed  at  200  GPa,  while 
transverse  modulus  varied  from  10  GPa  to  40  GPa,  and  (2)  when  the  transverse  modulus  was 
fixed  at  20  GPa,  while  longitudinal  modulus  varied  from  20  GPa  to  1000  GPa. 

Figure  67  shows  that  the  effective  moduli  of  the  carbon  foam  highly  depended  on 
the  change  of  transverse  modulus  of  the  foam  ligaments.  The  moduli  prediction  with  the 
anisotropic  ligament  properties  of  E^  =  200  GPa  and  E^  =  20  GPa  yielded  similar  bulk  tensile 
moduli  of  the  carbon  foam  with  the  isotropic  ligament  property  of  E-25  GPa  .  Meanwhile, 
Figure  68  shows  that  the  change  of  the  longitudinal  bulk  moduli  affected  the  foam  moduli  less 
significantly  than  the  change  of  the  transverse  ones  as  observed  in  Figure  67.  As  Figure  69 
shows,  the  rate  of  increment  of  the  bulk  foam  moduli  decreased  as  the  longitudinal  moduli  of  the 
ligament  increased  at  all  densities.  Higher  dependency  on  the  transverse  moduli  of  the  foam 
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Figure  67.  Prediction  of  Bulk  Tensile  Moduli  of  the  Carbon  Foam  Predicted  with  Isotropic  and 

Anisotropic  Properties  of  Foam  Ligaments.  Longitudinal  modulus  of  the  ligament  was 
fixed  at  200  GPa,  while  transverse  modulus  varied  from  10  GPa  to  40  GPa. 
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Figure  68.  Prediction  of  Bulk  Tensile  Moduli  of  the  Carbon  Foam  Predicted  with  Isotropic  and 

Anisotropic  Properties  of  Foam  Ligaments.  Transverse  modulus  of  the  ligament  was 
fixed  at  20  GPa,  while  longitudinal  modulus  varied  from  20  GPa  to  1000  GPa. 
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Figure  69.  Variation  of  Buik  Tensile  Moduli  of  the  Carbon  Foam  with  the  Increase  of  Longitudinal 
Modulus  of  Foam  Ligament.  Rate  of  increment  in  the  foam  moduli  decreased  as  the 
longitudinal  moduli  increased  at  all  densities. 
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Figure  70.  Prediction  of  Buik  Poisson’s  Ratio  of  the  Carbon  Foam  Predicted  with  isotropic  and 

Anisotropic  Properties  of  Foam  Ligaments.  Longitudinai  moduius  of  the  iigament  was 
fixed  at  200  GPa,  whiie  transverse  moduius  varied  from  10  GPa  to  40  GPa. 
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Figure  71.  Prediction  of  Bulk  Poisson’s  Ratio  of  the  Carbon  Foam  Predicted  with  Isotropic  and 
Anisotropic  Properties  of  Foam  Ligaments.  Transverse  modulus  of  the  ligament  was 
fixed  at  20  GPa,  while  longitudinal  modulus  varied  from  60  GPa  to  1000  GPa. 


ligaments  than  on  the  longitudinal  ones  illustrated  that  the  foam  ligaments  were  under  not  only 
axial  deformation  but  also  significant  bending  deformation.  Furthermore,  the  bending 
deformation  occurred  as  a  combination  of  pure  bending  and  transverse  shear  deformation,  which 
can  be  explained  with  Timoshenko’s  beam  theory. 

Figures  70  and  71  show  the  predicted  Poisson’s  ratio  of  the  carbon  foam  at 
various  densities.  As  in  the  case  of  the  isotropic  ligament  properties,  the  present  prediction 
yielded  nearly  constant  Poisson’s  ratios  unless  the  carbon  foam  was  of  low  density.  Unlike  the 
bulk  tensile  modulus,  the  Poisson’s  ratio  was  dependent  on  both  the  longitudinal  and  transverse 
moduli  of  the  foam  ligament,  as  shown  in  both  figures. 


7.3  SUMMARY  AND  RECOMMENDATIONS 

The  emerging  ultralightweight  material,  carbon  foam,  was  modeled  with  3-D 
microstructures  to  develop  a  basic  understanding  of  the  performance  of  open-cell  foam  materials. 
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To  develop  a  reliable  model  for  the  carbon  foam,  a  quantitative  microstructural  characterization 
of  foam  ligaments  and  nodes  should  be  complemented  with  the  model  development. 

Because  of  the  randomness  and  complexity  of  the  microstructure  of  the  carbon  foam,  the 
representative  cell  ligaments  were  first  characterized  in  detail  at  the  microstructural  level.  The 
microstructural  characteristics  (or  properties)  were  then  correlated  with  the  bulk  properties 
through  the  present  model.  A  series  of  databases  should  be  collected  for  various  size  and  spatial 
orientation  of  the  cell  ligaments,  as  well  as  the  property  variation  due  to  the  graphitic  alignment 
along  the  longitudinal  direction  of  the  ligaments.  The  model  is  thus  expected  to  provide  a  basis 
for  establishing  a  process-property  relationship  and  optimizing  foam  properties. 

The  present  model  yielded  a  fairly  reasonable  prediction  of  the  effective  bulk  properties 
of  the  foams.  The  present  predictions  were  compared  with  the  existing  simpler  model  and/or 
experimental  measurement.  We  observed  that  the  elastic  bulk  properties  of  the  foams  were 
dominated  by  the  bending  behavior  at  low- density  range  but  became  influenced  more  and  more 
by  the  axial  compliance  as  the  density  increased.  The  prediction  indicated  that  the  effective  bulk 
moduli  of  the  carbon  foam  were  dependent  more  on  the  transverse  modulus  of  the  foam 
ligaments  than  on  the  longitudinal  one.  The  rate  of  increment  in  the  foam  moduli  decreased  as 
the  longitudinal  moduli  increased  at  all  densities.  Unlike  the  Young’s  modulus,  the  Poisson’s 
ratio  was  dependent  on  both  the  longitudinal  and  transverse  moduli. 
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8.  SCANNING  ELECTRON  MICROSCOPY  INVESTIGATION  OE  DEFORMATION 
IN  THE  VICINITY  OF  A  MODE  II  CRACK  TIP  IN  UNIDIRECTIONAL 
TOUGHENED  EPOXY  COMPOSITES 

Delamination- related  research  on  fiber- reinforced  composites  has  attracted  considerable 
interest  due  to  its  significant  practical  importance.  Significant  differences  in  deformation 
patterns  in  the  crack  tip  vicinity  were  observed  for  Mode  I  and  Mode  II  delaminations  and 
reported  by  several  authors,  including  [45-47].  A  sequence  of  45°  (with  respect  to  the  crack 
plane)  microcracks  is  believed  to  form  beyond  the  Mode  II  crack  tip  in  polymer  composites.  The 
propagation  of  a  Mode  II  crack  is  believed  to  be  a  result  of  rupture  and  coalescence  of  these 
tension  microcracks. 

8.1  SPECIMEN  PREPARATION 

Two  types  of  specimens  were  utilized  for  testing.  A  24-ply  unidirectional  miniature 
specimen  (1.4”  long  and  0.25”  wide  )  was  cut  out  of  a  standard  size  (8”  long,  1”  wide)  specimen 
after  a  shear  crack  formed  in  the  course  of  end  notch  flexure  (ENF)  testing.  The  miniature 
specimen  was  cut  out  so  that  the  end  of  the  shear  crack  observed  in  the  optical  microscope  would 
be  halfway  between  the  bottom  support  and  the  bending  load  application  point.  The  material 
system  utilized  in  the  experiment  was  a  toughened  carbon  epoxy  unidirectional  composite, 
IM7/977-3. 

The  second  type  of  specimen  utilized  was  a  sandwich  plate  with  eight-ply  IM7/5250-4 
face  sheets  and  a  0.017”-thick  828  neat  epoxy  layer  in  between.  The  sandwich  plate  was 
manufactured  by  using  0.005”  spacers  on  each  sheet  filled  with  resin  and  then  assembled  into  the 
sandwich  at  room  temperature  cure.  A  Teflon  tape  was  inserted  between  the  two  assemblies  to 
create  a  starter  crack.  The  same  size  specimens  (1.4”  by  0.25”)  ^\ere  cut  out  from  the  sandwich 
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panel  so  the  Teflon  tape  would  end  halfway  between  the  bottom  support  and  the  bending  load 
application  point. 

8.2  EXPERIMENTAL  PROCEDURE 

A  Philips  XL30  environmental  scanning  electron  microscope  (E-SEM)  was  used  for 
testing.  A  three-point  bending  fixture  with  a  1000- lb  tensile  substage  for  XL30  E-SEM  was  used 
to  perform  in  situ  ENF  tests  on  miniature  specimens  described  in  the  previous  section.  The 
distance  between  the  lateral  supports  was  0.8”.  The  load  cell  was  equipped  with  force  and  head 
displacement  controller  with  digital  readout.  Compressive  load  was  applied  to  create  three-point 
transverse  bending  of  the  specimens. 

8.3  RESULTS  AND  DISCUSSION 

First,  the  24-ply  precracked  IM7/977-3  specimen  was  tested.  After  setting  up  the 
specimen  in  the  SEM,  it  was  loaded  in  bending  until  crack  extension  took  place.  The  examination 
of  the  crack  tip  zone  revealed  details  of  deformation  patterns,  which  at  first  appeared  similar  to 
those  reported  in  the  literature  [45,47].  A  typical  sequence  of  45 -degree  microcracks  was 
observed  in  the  crack  tip  area,  as  shown  in  Figure  72(a).  However,  as  can  be  seen  in  the  same 
figure,  the  damage  has  already  extended  beyond  this  microcracking  area.  It  appears  that  the 
microcracking  pattern  in  this  case  may  be  a  post- fracture  effect  when  a  thin  surface  layer  of 
polymer  remains  intact  after  delamination  and  then  ruptures  when  the  opening  displacement 
increases.  To  confirm  that  the  shadowed  area  on  the  right  side  of  the  figure  is  a  delamination, 
another  micrograph  was  taken  at  a  slightly  higher  load  and  shown  in  Figure  72(b).  A  clear 
delamination  can  be  seen  on  the  right  side  of  the  figure,  confirming  that  the  area  ahead  of  the 
microcracking  in  Figure  72(a)  was  a  delamination.  The  surface  layer  microcracks  ruptured  and 
coalesced  but  after  the  Mode  II  delamination  had  already  propagated. 
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(a) 


(b) 

Figure  72.  SEM  Image  of  the  Mode  II  Crack  Tip  Region  at  Two  Consecutive  Load  Levels. 


The  behavior  observed  in  the  present  work  can  be  attributed  to  the  following  factors.  The 
resin  system  utilized  in  the  experiments  is  a  toughened  epoxy,  which  can  exhibit  very  ductile 
behavior  and  form  surface  films  behaving  as  described  above.  However  it  is  believed  that  the 
main  difference  between  the  composite  tested  in  the  present  work  and  that  reported  in  [45-47]  is 
the  absence  of  the  resin-rich  area  where  the  Mode  II  crack  propagated  in  the  present  study. 
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A  sandwich-type  specimen  containing  Teflon  tape  to  start  the  interlaminar  crack  in  the 
resin-rich  region  was  tested  next.  Figure  73  shows  the  image  of  the  region  inside  the  resin  layer 
near  the  end  of  the  Teflon  tape  under  135  lbs  of  transverse  load  with  the  starter  crack  just 
emanating  from  the  edge  of  the  Teflon  sublayer.  One  can  observe  several  fragments  of  the  crack 
parallel  to  the  face  sheet  plane.  It  was  determined  that  the  plane  of  crack  propagation  coincided 
with  the  artificial  interface  created  in  the  center  of  the  adhesive  layer  during  manufacture  of  the 
sandwich.  This  interface  was  created  by  assembling  the  two  face  sheets  with  the  828  resin  spread 
on  each  face  sheet.  The  shear  crack  continued  to  extend  under  loading  up  to  186  lb  until 
deflecting  from  the  plane  toward  the  outer  layer  and  forming  a  cavity,  shown  in  Figure  74(a).  The 
cavity  continued  to  increase  [Figure  74(b)];  however,  no  longitudinal  extension  of  the  crack  took 
place.  This  behavior  is  consistent  with  the  fact  that  the  transverse  shear  stress  vanishes  at  the 
center  (under  the  bending  fixture)  of  the  specimen  at  its  midplane.  It  is  also  consistent  with  the 
crack  termination  location  from  specimen  to  specimen;  a  total  of  three  replicas  were  investigated. 


Aoc  V  Spol  Magn  Del  WO  I - 1  200  |im 

IS  O  kV  6.0  I26x  GSEO.e  S.6  Toti  82S/D230  RESIN  LAVER  06280 


Figure  73.  SEM  Image  of  the  Mode  II  Crack  Tip  Region  in  the  Resin-Rich  Region  of  the  Sandwich 
Specimen. 
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Figure  74.  Fully  Extended  Interlaminar  Crack  at  (a)  186-lb  Loading  and  (b)  226-lb  loading. 
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Further  loading  led  to  fiber  failure,  including  compression  failure  near  the  inner  surface  of  the 
lower  face  sheet,  as  shown  in  Figure  75.  It  is  interesting  to  note  that  only  distributed  fiber 
breakage  was  observed  on  the  surface  of  the  lower  face  sheet,  which  was  loaded  in  tension. 

8.4  SUMMARY  AND  RECOMMENDATIONS 

In  situ  observation  of  deformation  and  failure  phenomena  in  advanced  composite 
materials  was  conducted  in  a  high-resolution  Philips  XL30  E-SEM.  Interlaminar  crack  tip 
features  were  examined  in  unidirectional  fiber- reinforced  composites  with  and  without  the  resin- 
rich  area  at  the  midsurface  of  the  laminate.  Face  sheet  failure  was  observed  in  the  sandwich 
plates  at  relatively  small  extension  of  the  Mode  II  crack  in  the  adhesive.  Further  work  is  required 
to  investigate  the  Mode  II  crack  extension  mechanism  in  a  resin-rich  region. 


Figure  75.  Compression  Fiber  Faiiure  on  Inner  Surface  of  the  Bottom  Face  Sheet  at  226-lb  Load. 
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9.  NONLOCAL  MODELS  OF  STRESS  FIELD  CONCENTRATIONS  AND  EFFECTIVE 
THERMOELASTIC  PROPERTIES  OF  RANDOM 
STRUCTURE  COMPOSITES  [48] 

In  this  section  we  consider  a  linearly  thermoelastic  composite  medium,  which  consists  of  a 
homogeneous  matrix  containing  a  statistically  homogeneous  set  of  ellipsoidal  inclusions.  The 
multiparticle  effective  field  method  (MEFM)  (see  for  reference  [49])  based  on  the  theory  of 
functions  of  random  variables  and  Green's  functions  is  used.  Within  this  method  a  hierarchy  of 
statistical  moment  equations  for  conditional  averages  of  the  stresses  in  the  inclusions  is  derived. 
The  hierarchy  is  established  by  introducing  the  notion  of  an  effective  field.  In  this  way  the 
interaction  of  different  inclusions  is  taken  directly  into  account  in  the  framework  of  the 
homogeneity  hypothesis  of  the  effective  field.  Combining  the  MEFM  with  the  standard  scheme 
of  the  iteration  method  (IM)  and  Fourier  transform  method  (FTM)  permits  one  to  obtain  the 
explicit  representations  for  the  nonlocal  integral  and  differential  operators,  respectively,  of  any 
order  describing  overall  effective  properties  as  well  as  the  stress  concentrator  factor  in  the 
components.  It  is  shown  that  the  reduction  of  the  integral  operator  to  the  differential  one  for 
sufficiently  smooth  statistical  average  stress  fields  is  superior  to  the  FTM.  After  some  additional 
assumptions  the  method  proposed  is  reduced  to  the  perturbation  method  as  well  as  to  the 
quasicrystalline  approach.  For  some  concrete  numerical  examples,  one  can  demonstrate  the 
advantage  of  the  IM  over  the  FTM. 

The  prediction  of  the  behavior  of  composite  materials  by  the  use  of  mechanical  properties 
of  constituents  and  their  microstructure  is  a  central  problem  of  micromechanics,  which  is 
eventually  reduced  to  the  estimation  of  stress  fields  in  constituents.  A  considerable  number  of 
methods  are  known  in  the  linear  theory  of  statistically  homogeneous  composites  being  considered 
in  this  paper  which  yield  the  effective  elastic  constants  and  stress  field  averages  in  the 
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components.  Appropriate,  but  by  no  means  exhaustive,  references  are  provided  by  the  reviews 
[48,50] .  When  the  statistical  average  stress  in  a  composite  medium  varies  sufficiently  slowly 
compared  to  the  size  scale  of  the  microstructure  (i.e.,  roughly  speaking,  the  separation  of  the 
external  and  internal  scales  takes  place),  the  material  macroscopically  behaves  as  a  homogeneous 
body  with  some  effective  moduli,  which  can  be  estimated  by  the  different  methods  referred  to 
above.  If  the  condition  of  the  separation  of  scales  summarized  above  is  not  valid,  the  material's 
response  cannot  be  adequately  described  in  the  framework  of  the  local  theory  of  elasticity  for  the 
homogeneous  medium. 

For  the  analysis  of  the  widely  separated  scales  (but  not  too  widely),  a  number  of 
phenomenological  approaches  have  been  proposed  to  enhance  the  continuum  model  by  nonlocal 
terms,  either  introduced  through  an  integral  equation  or  through  an  additional  gradient  equation 
based  on  the  assumption  that  the  forces  between  material  points  can  be  long  range  in  character 
(see  for  reference  [48]).  The  different  modifications  of  these  phenomenological  constitutive 
relations  were  used  for  the  analyses  of  various  size  effects,  such  as  strain  gradient  hardening,  size 
effect  in  torsion  and  indentation,  strain  gradient  plasticity,  and  nonlocal  damage  mechanics. 

The  micromechanical  models  are  used  to  answer  a  fundamental  question  of  how  length 
scales  in  the  effective  constitutive  equations  could  be  directly  derived  from  the  geometrical  and 
mechanical  properties  of  constituent  phases.  It  is  known  that  the  eventual  abandonment  of  this 
hypothesis  of  statistically  homogeneous  fields  leads  to  a  nonlocal  coupling  between  statistical 
averages  of  the  stress  and  strain  tensors  when  the  statistical  average  stress  is  given  by  an  integral 
of  the  field  quantity  weighted  by  some  tensorial  kernel,  i.e.,  the  nonlocal  effective  elastic 
operator.  Analysis  of  highly  contrasted  statistically  homogeneous  media  is  simplified  for 
sufficiently  smooth  external  loading  permitting  the  use  of  the  Taylor  expansion  for  the  statistical 
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average  of  a  stress  field  and  the  application  of  an  FTM.  In  so  doing  the  initial  integral  equation  is 
reduced  to  an  algebraic  polynomial  equation  with  constant  coefficients  in  a  Fourier  space  with 
forthcoming  solution  and  the  implementation  of  the  inverse  Fourier  transform.  The  scheme 
summarized  informally  above  is  usually  based  on  the  hypothesis  of  the  homogeneity  of  the 
effective  field  in  which  each  particle  is  located.  The  quasicrystalline  approximation  by  Lax  is 
often  used  for  truncation  of  the  hierarchy  of  integral  equations  involved,  leading  to  neglect  of 
direct  multiparticle  interactions  of  inclusions.  This  reduced  the  analysis  to  the  use  of  statistical 
information  of  up  to  two-point  correlation  functions  and  allowed  one  to  derive  explicit  relations 
for  the  nonlocal  overall  differential  operator  by  the  different  methods:  by  the  effective  field 
method  or  by  the  method  of  conditional  moments  (see  for  reference  [48]). 

The  main  disadvantage  of  the  quasicrystalline  approximation,  consisting  of  the  neglect  of 
direct  multiparticle  interactions  of  inclusions,  was  overcome  recently  by  the  MEFM  (references 
may  be  found  in  the  survey  [49]).  The  MEFM  is  based  on  the  theory  of  functions  of  random 
variables  and  Green's  functions.  Within  this  method  a  hierarchy  of  statistical  moment  equations 
for  conditional  averages  of  the  stresses  in  the  inclusions  is  derived.  The  hierarchy  is  established 
by  introducing  the  notion  of  an  effective  field.  In  this  way  the  interaction  of  different  inclusions 
is  taken  directly  into  account  in  the  framework  of  the  homogeneity  hypothesis  of  the  effective 
field.  Combining  the  MEFM  with  the  standard  scheme  of  the  FTM  usually  used  for  the  analysis 
of  the  statistically  inhomogeneous  stress  fields  allowed  (see  [51,52])  to  obtain  the  explicit 
representation  of  a  nonlocal  overall  operator  in  the  form  of  the  second -order  differential  operator. 
The  author  ]53]  estimated  the  nonlocal  integral  operator  of  the  overall  constitutive  equation  for 
periodic  structure  composites  by  the  use  of  the  first  iteration  of  the  IM  and  showed  the  reduction 
of  the  integral  operator  obtained  to  the  differential  operator  involved  derived  by  the  FTM. 
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The  last  approach  is  generalized  in  this  report  for  the  analysis  of  statistically  homogeneous 
random  structure  composites.  By  the  use  of  the  IM  and  FTM,  one  obtains  the  explicit 
representations  for  the  nonlocal  integral  and  differential  operators,  respectively,  of  any  order 
describing  overall  effective  properties  as  well  as  the  stress  concentrator  factor  in  the  components. 
We  show  the  reduction  of  the  integral  operator  to  the  differential  one  for  sufficiently  smooth 
statistical  average  stress  fields.  After  some  additional  assumptions  the  method  proposed  is 
reduced  to  the  perturbation  method  as  well  as  to  the  quasicrystalline  approach.  In  some  concrete 
numerical  examples,  one  demonstrates  the  advantage  of  the  IM  over  the  FTM. 

9.1  THE  GENERAL  SCHEME 

Let  a  linear  elastic  body  occupy  an  open-bounded  domain  with  a  smooth  boundary  and 
contain  a  homogeneous  matrix  and  a  statistically  homogeneous  set  X  of  ellipsoidal  inclusions  v,. 

It  is  assumed  that  the  inclusions  can  be  grouped  into  components  (phases)  with  identical 
mechanical  and  geometrical  properties  (such  as  the  shape,  size,  orientation,  and  microstructure  of 
inclusions).  For  the  sake  of  definiteness,  in  the  2-D  case  we  will  consider  a  plane-strain  problem. 
At  first  no  restrictions  are  imposed  on  the  elastic  moduli  L(x)  of  phases. 

To  simplify  the  exact  infinite  system  of  integral  equations  obtained,  we  now  apply  the 
main  hypothesis  of  many  micromechanical  methods,  the  approximation  called  the  effective  field 
hypothesis: 

HI)  Each  inclusion  v,  has  an  ellipsoidal  form  and  is  located  in  the  homogeneous  effective 
field  which  is  found  from  self-consistent  estimations. 

H2)  For  a  sufficiently  large  number  of  interacting  inclusions  n,  the  effective  fields  acting 
on  the  n  and  n+1  inclusions  coincide. 
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The  fundamental  differences  of  hypothesis  H2  and  the  quasicrystalline  approximation  by  Lax 
were  discussed.  According  to  hypothesis  HI  and  in  view  of  the  linearity  of  the  problem,  there 
exist  constant  fourth-  and  second-rank  tensors  describing  the  average  response  of  n  inclusions 
subjected  to  the  effective  field. 

Then  the  infinite  system  of  integral  equations  can  be  reduced  to  single  integral  equations 

(45) 

for  statistical  average  of  stresses  in  the  inclusions  <o>;,  where  Y  is  a  tensor  and  K  is  an  integral 
operator. 

Thus  we  reduced  the  infinite  system  of  integral  equations  to  the  standard  form  of  operator 
Eq.  (45)  with  the  regular  integral  kernel  of  convolution  type  that  can  be  solved  by  such  known 
methods  as  the  method  of  mechanical  quadratures,  successive  approximations,  and  FTM. 
Application  of  the  method  of  mechanical  quadratures  is  very  popular,  although  it  is  required  to 
solve  a  linear  system  of  very  high  order  even  for  a  smoothly  varying  load.  Because  of  this,  at  first 
we  will  solve  Eq.  (45)  by  the  method  of  successive  approximations,  which  is  also  called  the 
Neumann  series  method.  Their  connection  with  the  FTM  is  considered  in  [48]. 

The  IM  is  based  on  the  recursion  formula, 

,  (46) 

to  construct  a  sequence  of  functions  that  can  be  treated  as  an  approximation  of  the  solution 

of  Eq.  (45).  In  effect  the  IM  (Eq.  46)  transforms  the  integral  equation  problem  (Eq.  45)  into  the 
linear  algebra  problem  in  any  case.  Usually  the  driving  term  of  this  equation  is  used  as  an  initial 
approximation. 

Now  we  will  consider  the  FTM.  One  assumes  that  the  statistical  average  of  stress 
polarization  tensor  belongs  to  the  class  of  m  times  continuously  differential  functions.  Since  we 
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desire  an  explicit  representation  for  the  stress  concentrator  factor,  we  will  approximate  stresses  by 
the  first  terms  of  its  Taylor  expansion  and  reduce  Eq.  (45)  to  the  following  symbolic  equation: 

P (x,?)  <o>i  (x).  (47) 

The  operator  P(x,?)  is  a  linear  partial  differential  operator  with  the  constant  coefficients 
for  statistically  homogeneous  media.  Considering  that  Eq.  (47)  is  a  differential  equation  with 
constant  coefficients,  the  method  of  solution  that  first  comes  to  mind  is  using  the  ETM  to 
transform  the  differential  problem  of  solving  Eq.  (47)  into  the  division  problem  of  solving  the 
multiplicative  equation.  Then  the  solution  of  Eq.  (47)  is  represented  by  the  differential  operator 

<a>,-  (x)=Q(?)a®''®''®s®‘^  (x),  (48) 

where  an  explicit  form  of  the  differential  operator  Q(?)  of  an  infinite  order  is  presented  in  [48] . 

We  obtained  the  representations  of  integral  and  differential  operators  of  an  infinite  order 
for  the  stress  concentrator  factors  as  well  as  for  the  effective  thermoelastic  properties.  The 
relations  obtained  were  simplified  for  some  particular  cases  such  as  weakly  inhomogeneous 
media  as  well  as  the  quasicrystalline  approximation.  We  demonstrated  the  reduction  of  integral 
operators  to  the  differential  ones. 

9.2  NUMERICAL  RESULTS 

Just  for  demonstration  of  the  comparison  of  available  experimental  data  with  the 
predicting  opportunity  of  the  method  proposed,  we  will  consider  the  zero- order  approximation  of 
the  method  that  is  the  estimation  of  the  effective  elastic  moduli  L*.  Assume  the  matrix  is  epoxy 
resin  (^*^^^=4.27  GPa  and  p*^°^=1.53  GPa)  which  contains  identical  circular  glass  fibers  (k^^^=50.S9 
GPa  and  fJ ^^=35. 04  GPa).  Two  alternative  radial  distribution  functions  g( r)  of  inclusion 
distribution  will  be  examined:  the  step  function  and  the  function  proposed  in  [54]  that  takes  into 
account  a  neighboring  order  in  the  distribution  of  the  inclusions.  As  can  be  seen  in  Figure  76,  the 
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Figure  76.  Variation  of  the  Effective  Shear  Moduius  |x*  as  a  Function  of  a  Concentration  of  the 

inclusions  c.  Experimental  data  (circles)  and  curves  calculated  by  MEFM  and  nonstep 
function  g(r)  (solid  line),  MEFM  and  step  function  g(r)  (dot-dashed  line),  perturbation 
method  and  nonstep  g(r)  (dashed  line),  perturbation  method  and  step  function  g(r) 
(diamonds),  Mori-Tanaka  approach  (dotted  line). 

use  of  the  approach  based  on  the  quasicrystalline  approximation  [also  called  Mori-Tanaka  (MT) 
approach]  leads  to  an  underestimation  of  the  effective  shear  modulus  by  1.85  times  for  c=0.7 
compared  to  the  experimental  data  as  well  as  to  the  more  exact  approximation  of  the  MEFM, 
which  provides  the  best  comparison  with  experimental  data  (see  [55]).  It  should  be  mentioned 
that  the  estimations  by  the  MEFM  are  slightly  sensitive  to  the  choice  of  the  radial  distribution 
function  g.  It  was  observed  that  the  proposed  formulation  compared  well  with  experimental  data 
for  c  up  to  65  percent;  for  fiber  volume  fractions  greater  than  70  percent,  microdefects  existing  in 
experimental  specimens  may  significantly  affect  the  overall  elastic  moduli.  At  the  same  time,  the 
MT  solution  differs  from  the  perturbation  method  approximation  by  not  more  than  five  percent 
for  the  concentration  of  the  inclusions  c<0.7  and  does  not  depend  on  the  radial  distribution 
function  (in  contrast  to  the  perturbation  method). 

Our  approach  to  addressing  this  section  will  be  to  employ  the  nonlocal  equations  for  stress 
concentrator  tensors  we  derived  in  the  previous  sections  by  the  IM  and  FTM.  We  will  consider 
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ensemble -averaged  stresses,  and  determine  the  stress  concentrator  factors  at  the  inclusions  by  the 
IM  and,  in  some  particular  cases,  by  the  FTM.  The  quantitative  results  we  obtain  will  be  for  the 
2-D  case  (plane  strain  problem)  for  the  two -phase  composites  consisting  of  an  isotropic  matrix 
reinforced  by  a  random  dispersion  of  isotropic  identical  circle  particles. 

Let  us  now  demonstrate  the  application  of  the  theoretical  results  by  considering  an 
isotropic  composite  made  of  an  incompressible  isotropic  matrix  filled  with  rigid  circle  inclusions 
of  one  size.  This  example  was  chosen  deliberately  because  it  provides  the  maximum  difference 
of  predictions  of  effective  elastic  response  estimated  by  various  methods,  and  was  considered  by  a 
number  of  authors.  We  shall  consider  the  response  for  a  normal  loading  that  varies  with  position 
in  its  loading  direction  <s  ii>=f(xi)diidijWi\h  three  specific  cases  of  the  functions 


fi(xi)=sin(pxi/4), 

(49) 

f2(x  1  )=0.6579\x  1 P  X 

(50) 

f3(X])=0.6584\x]\^-^^‘^  e  x/, 

(51) 

f4(xi)=0.6580xi^e 

(52) 

The  loading  with/i  belonging  to  (R)  was  analyzed  by  the  FTM  in  detail  for  the  different 
arguments  and  the  different  concentrations  of  the  spherical  inclusions.  We  can  only  define  the  m- 
th  derivatives  of  the  functions ^(0)  (j=2,3)  if/'"Y-0)=/'"V +0).  Now/2  belongs  to  and/j 
belongs  to  C^,  and  the  third  and  the  second  derivatives  do  not  exist  for  the  functions /2  and/j, 
respectively.  The  functions/  (j=2,3,4)  have  approximately  the  same  maximum,  and  the  same 
maximum  of  their  first  derivatives,  as  the  function  fi  (see  Figure  77).  In  so  doing,  independently 
from  the  concrete  micromechanical  average  scheme,  the  FTM  will  predict  zero  nonlocal  effects  at 
the  point  0  for  the  function /2(xi).  Analogous  analysis  leads  to  a  more  dramatic  conclusion  for  the 
function 
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Figure  77.  The  Functions  and  their  First  Derivatives  versus  Argument  f2(Xi)  (solid  line),  fi(xi) 
(dot-dashed  line),  f/Cxr.)  (dashed  line),  (dotted  line). 

fs,  that  is,  the  FTM  cannot  be  applied  in  principle  to  the  field  described  by  fs  because  the 
differentiable  function/jfx;  j,  which  is  very  close  to  f2(xi),  is  not  twice  differentiable  at  xi=0. 

The  comparative  results  estimated  by  the  MEFM  for  the  nonstep  radial  distribution 
function  and  for  different  functions (i=l,2,3,4)  for  the  0th  and  7th  iterations  are  presented  in 
Reference  [48]',  the  results  for  functions/;  differ  from  one  another  by  less  than  one  percent. 

The  first  few  iterations  of  estimations  of  stresses  for  the  function  /2  are  presented  in  Figure 
78.  As  can  be  seen,  a  few  order  approximations  of  the  IM  lead  to  significantly  different 
normalized  stresses  provided  by  the  derivatives/'”^  (if  they  exist),  whereas  the  second-order 
approximation  of  the  FTM  leads  to  the  degenerate  results. 

hi  [48]  only  the  function /2  is  analyzed  in  the  framework  of  the  iteration  scheme.  The 
comparative  analysis  of  the  MEFM,  MT  method,  and  the  perturbation  method  for  the  different 
radial  correlation  functions  is  presented.  As  is  shown,  the  MEFM  is  most  sensitive  to  the  choice 
of  radial  distribution  function  and  leads  to  the  maximum  nonlocal  effects  predicted.  It  is 
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Normalized  coordinate 

Figure  78.  The  First  Few  iterations  of  Statisticai  Averages  <Srr>;Y>frj versus  Xj  Estimated  for  the 
Functions  f2(Xi)  and  Nonstop  g(r)  by  the  MEFM.  Zero-order  (soiid  iine),  first-order 
(dotted  line),  second-order  (dot-dashed  line),  seventh-order  (dashed  line) 
approximations. 

interesting  that  nonlocal  response,  in  contrast  to  the  local  one  estimated  by  the  MT  method, 
depends  on  the  radial  distribution  function  g( r). 

9.3  SUMMARY  AND  RECOMMENDATIONS 

Let  us  discuss  the  main  scheme  as  well  as  the  short  sketch  of  limitations  and  of  possible 
generalization  and  application  of  the  methods  proposed. 

The  obtained  relations  depend  on  the  values  associated  with  the  mean  distance  between 
inclusions,  and  do  not  depend  on  the  other  characteristic  size,  i.e.,  the  mean  inclusion  diameter. 
This  fact  may  be  explained  by  the  initial  use  of  hypothesis  HI  dealing  with  the  homogeneity  of 
the  effect  inside  each  inclusion.  In  the  case  of  a  variable  representation  of  the  effective  field,  for 
instance,  in  polynomial  form,  the  mean  size  of  the  inclusions  will  be  contained  in  the  nonlocal 
dependence  of  microstresses  on  the  average  stress.  Such  an  improvement  based  on  the 
abandonment  from  the  hypothesis  HI  was  schematically  considered  in  [49]. 
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It  should  be  mentioned  that  the  IM  and  the  FTM  have  a  series  of  advantages  and 
disadvantages,  and  it  is  crucial  for  the  analyst  to  be  aware  of  their  range  of  application.  The  IM 
has  two  known  drawbacks.  First,  the  Neumann  series  ensures  the  existence  of  solutions  to 
integral  equations  of  the  second  kind  only  for  sufficiently  small  kernels  involved,  and  second,  in 
general,  it  cannot  be  summed  in  closed  form.  Of  course  Eq.  (45)  can  be  solved  directly  by  the 
quadrature  method  even  if  the  condition  of  smallness  of  the  kernels  is  not  valid.  However, 
strongly  inhomogeneous  problems  may  lead  to  much  larger  numbers  of  quadrature  points,  making 
iteration  potentially  worthwhile.  Moreover,  increasing  the  problem  dimensionality  (from  2-D  to 
3-D)  raises  the  number  of  nodes  to  the  dimensional  power,  and  the  situation  changes  radically. 

As  was  shown,  only  a  few  iterations  are  necessary;  these  iterations  prove  very  much  faster  than  a 
direct  inversion  of  the  nonlocal  operator  by  the  quadrature  method. 

The  reduction  of  integral  operators  to  the  differential  ones  allows  an  understanding  of 
drawbacks  of  the  FTM.  The  first  one  is  that  for  obtaining  an  m-order  differential  operator,  it  is 
necessary  that  s  belongs  to  C'".  In  so  doing,  the  IM  providing  the  accuracy  of  differential 
operator  of  an  infinite  order  does  not  have  even  continuity  of  s  since  integration  is  a 

smoothing  operation  and  the  right-hand- side  integral  (Eq.  45)  is  likely  to  be  a  rather  smooth 
function  even  when  s  is  very  jagged.  However,  even  the  m-th  approximation  of  the  IM 
contains  the  differential  operators  of  infinite  orders  that  are  lost  in  the  FTM.  The  question  of  the 
convenience  of  using  one  method  over  another  is  solved  also  in  favor  of  the  IM,  because  in  the 
FTM  it  is  necessary  to  calculate  the  cumbersome  tensors  and  completeness  of  estimation,  which 
increases  dramatically  with  m,  while  in  the  IM  it  is  enough  to  estimate  a  single  tensor  and 
consecutively  to  apply  the  recursion  scheme  (Eq.  46),  the  completeness  of  which  does  not  depend 
on  the  iteration  number.  Thus,  the  single  advantage  of  the  FTM  comprised  of  obtaining 
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analytically  explicit  relations  dwindles  in  light  of  the  disadvantages  mentioned  above,  such  as  a 
requirement  of  smoothness  of  s  belongs  to  C"*,  and  the  intricacy  of  analytical  calculations. 
We  indicated  only  mathematical  and  computational  difficultie  s  in  the  use  of  the  FTM,  which  can 
be  solved,  at  least  in  principle,  at  a  cost  of  great  efforts  if  the  analytical  solution  is  necessary. 

However,  it  should  be  mentioned  that  there  is  an  extremely  important  class  of 
micromechanical  problems  for  statistically  inhomogeneous  media  (such  as  functionally  graded 
and  clustered  materials),  analysis  of  which  by  the  FTM  is  questionable  in  principle.  The 
breakdown  of  the  assumption  of  statistical  homogeneity  leads  to  the  inequality  Y(x)?const.  Then 
the  average  stresses  are  not  a  constant.  However,  even  this  simplest  case  of  homogeneous 
boundary  conditions  leads  to  a  fundamental  prohibitive  obstacle  against  the  use  of  the  FTM. 
Indeed,  the  inhomogeneity  of  average  stresses  yields  the  fact  that  the  linear  differential  operator 
(Eq.  47)  has  variable  coefficients.  Then  the  jump  from  Eq.  (47)  to  Eq.  (48),  respectively,  based 
on  the  properties  of  the  Fourier  transform,  is  difficult,  and  an  applicability  of  the  FTM  is 
questionable.  In  so  doing,  the  use  for  statistically  inhomogeneous  media  of  the  IM  inserts 
requires  slight  modification  of  the  scheme  presented  [Eqs.  (45)  and  (46)]  (see,  e.g.,  [49])  for 
application  of  the  IM  for  research  of  periodic  graded  composites).  However,  the  analysis  of  these 
problems  is  beyond  the  scope  of  the  current  study  and  will  be  considered  in  forthcoming 
publications  by  the  authors. 


121 


10.  MICROMECHANICS  OF  NONLOCAL  EFFECTS  IN  GRADED  MATERIALS 

A  considerable  number  of  methods  are  known  in  the  linear  theory  of  statistically 
homogeneous  eomposites;  appropriate,  but  by  no  means  exhaustive,  referenees  are  provided  by 
the  reviews  [49,50] .  Signifieantly  less  research  has  been  performed  for  statistically 
inhomogeneous  composites.  In  such  a  case,  the  ergodicity  fails,  and  ensemble  and  averages  do 
not  coincide.  The  degenerate  ease  of  this  material  is  a  random  matrix  composite  bonded  in  some 
directions  as  well  as  the  composite  medium  for  which  the  inclusions  are  located  in  a  region 
bonded  in  some  direetions,  although  unrestrietedness  of  the  domain  of  inclusion  locations  does 
not  preelude  statistieal  inhomogeneity.  For  example,  any  laminated  composite  material  randomly 
reinforced  by  aligned  fibers  in  each  ply  is  a  statistically  inhomogeneous  material.  The  authors 
[56]  proposed  statistieal  deseriptions  for  particulate,  statistically  inhomogeneous  two -phase 
random  media  by  the  use  of  the  theory  of  a  general  Poisson  process  and  estimated  for  some 
simulated  fully  penetrable  (Poisson  distribution)  spheres  the  eanonical  n-point  microstructure 
funetion,  the  nearest- neighbor  function,  and  the  linear-path  funetion  that,  unlike  the  homogeneous 
case,  will  depend  on  the  absolute  positions  of  their  argument. 

The  eoncept  of  elusters  is  similar  to  that  of  fraetal  struetures,  and  the  role  of  statistieal 
descriptor  can  be  treated  by  such  parameters  as  cluster  size,  fractal  dimension,  and  the  radius  of 
gyration.  Even  though  these  and  some  other  parameters  are  used  for  the  identification  of 
clustered  struetures,  it  is  not  suffieient  to  predict  the  overall  properties  of  composites  due  to 
several  reasons.  The  first  one  is  that  these  parameters  are  not  complete  enough  for  the 
eharaeterization  of  the  mieromorphology  of  fillers  simply  beeause  one  can  present  other 
morphology  with  the  same  descriptors.  More  informative  characteristics  of  the  random 
eonfigurations  use  statistical  second-order  quantities,  whieh  examine  the  assoeiation  fillers 
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relative  to  other  particles  in  an  immediate  local  neighborhood  of  the  reference  filler.  Another 
reason  is  that  the  prediction  of  mechanical  properties  requires  one  or  another  micromechanical 
model.  Advanced  micromechanical  models  normally  use  two  different  sorts  of  statistical 
descriptors  such  as  n-point  probability  density  and  an  alternative  descriptor  of  location  of  m 
inclusions  with  the  c&ai&[sxi,X2,...Xm  described  by  the  probability  density  (^m.(xi,X2,...Xm). 

The  interest  in  one  or  another  cluster  shape  is  defined  by  the  physical  problem  being 
considered.  If  the  problem  is  recrystallization  or  residual  stresses,  it  is  reasonable  that  the 
important  clusters  have  a  convex  shape.  If  one  is  concerned  by  the  damage  properties  of 
dielectric  breakdown,  a  long  chain  of  broken  or  conducting  particles  forming  a  percolating 
network  can  be  more  deleterious  than  the  same  number  of  particles  arranged  in  a  convex  cluster. 
For  chain  clusters,  not  only  the  local  density  of  particles  but  also  their  arrangement  is  relevant  that 
are  demonstrated  as  in  a  limiting  case  in  a  concept  of  weakest  link.  Due  to  the  numerous  possible 
cluster  configurations  and  limited  opportunity  of  prediction  models,  we  will  consider  the 
degenerate  case  of  cluster  materials  that,  nevertheless,  has  wide  applications.  Namely,  we  will 
analyze  so-called  ideal  cluster  materials  (see,  e.g.,  [57])  that  contain  either  finite  or  infinite, 
deterministic  or  random  ellipsoidal  domains  called  particle  clouds  distributed  in  the  composed 
matrix.  In  so  doing  the  concentration  of  particles  is  a  piecewise  and  homogeneous  one  within  the 
areas  of  ellipsoidal  clouds  and  composed  matrix.  In  particular,  in  this  section  we  consider  a 
single  particle  cloud  with  the  shape  of  a  thick  ply  located  in  an  infinite  matrix  with  zero 
concentration  of  particles.  The  ideal  cluster  configuration  may  be  created  by  distributing  cloud 
centers  deterministically  or  randomly  and  placing  offspring  points  randomly  and  uniformly 
around  parent  points  within  a  predetermined  area.  For  a  description  of  the  cluster  arrangement, 
the  notion  of  cluster  radius  is  used  for  the  presentation  of  a  cluster  as  a  particle  of  the  second 
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generation  with  the  radius  defined  by  the  average  distances  between  the  particles.  Usually  the 
cluster  covers  80  to  90  percent  of  the  total  volume  of  the  sphere  with  the  cluster  radius. 

In  many  cases,  the  spatial  arrangement  of  the  fibers  is  neither  uniformly  random  nor 
periodic.  For  a  statistical  description  of  the  microstructure  of  ceramic  matrix  composites 
containing  fiber- rich  and  fiber-poor  regions  that  have  a  laminated  structure,  the  authors  [58]  have 
proposed  the  so-called  strip  model.  The  key  parameters  of  the  strip  model  are  the  number  of 
densities  of  fiber  in  the  fiber-poor  and  fiber-rich  regions,  the  widths  of  some  fiber-poor  and  fiber- 
rich  regions,  and  the  global  number  of  densities  of  the  fibers  is  defined  by  the  a\€raging  of  these 
parameters.  Thus,  the  number  of  densities  in  a  strip  model  is  described  by  a  piecewise  constant 
function.  Generalization  of  this  model  is  based  on  the  introduction  of  a  variable  number  of 
densities  n(x2),  defined  as  the  average  number  of  disk  centers  per  unit  in  a  thin  ply  parallel  to  the 
axis  X2=0. 

In  this  section  we  consider  a  linearly  elastic  composite  medium,  which  consists  of  a 
homogeneous  matrix  containing  statistically  inhomogeneous,  so-called  graded  random  field  of 
inclusions.  For  functionally  graded  materials,  when  the  concentration  of  the  inclusions  is  a 
function  of  the  coordinates,  the  micromechanical  approach  is  based  on  the  generalization  of  the 
MEFM,  previously  proposed  for  statistically  homogeneous  random  structure  composites  by  the 
author  (see  for  reference  and  details  [49]).  The  nonlocal  integral  effective  operator  of  elastic 
effective  properties  is  estimated.  The  nonlocal  dependencies  of  the  effective  elastic  moduli,  as 
well  as  of  conditional  averages  of  the  strains  in  the  components,  on  the  concentration  of  the 
inclusions  in  a  certain  neighborhood  of  points  considered  are  detected. 
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10.1  GENERAL  SCHEME  OF  THE  APPROACH  PROPOSED 


The  trivial  generalization  of  the  approach  presented  by  Eq.  (45)  leads  to  the  estimation  of 
statistical  average  of  strains  inside  the  inclusions  as  well  as  the  effective  properties: 

<a>;(x)=Y(x)a‘‘''®’^s®(x)+K(x)<a>/ .  (53) 

The  relations  are  the  first-order  approximation  of  the  nonlocal  theory.  The  construction  of  the 
next  approximation  is  obvious.  Only  at  first  glance,  Eq.  (53)  is  equivalent  to  the  corresponding 
ones  obtained  for  global  effective  properties.  The  main  difference  is  that  Y(Ar)  and,  therefore, 
L*(x)  depend  on  the  parameters  of  the  inclusion  distribution  not  only  at  the  point  x,  but  also  in  a 
certain  neighborhood  of  that  point  leading  to  a  so-called  nonlocal  effect,  though,  of  course,  the 
effective  parameters  L*(x)  are  the  local  ones  in  the  sense  of  the  nonlocal  elasticity  theory.  The 
diameter  of  this  region  mentioned  above  is  estimated  as  three  times  the  characteristic  dimension 
of  the  inclusions.  As  a  result,  a  statistically  inhomogeneous  composite  medium  behaves  like  a 
macroscopically  inhomogeneous  medium  with  local  effective  modulus  L*(x)  determined  for  a 
nonlocal  distribution  of  the  inclusions  in  a  certain  neighborhood  of  the  point  considered. 

10.2  PREVIOUS  NUMERICAL  RESULTS 

Let  us  consider  a  strip  model  of  ideal  fiber  cluster  with  probability  densities  ^)i(x)=n  and 

(^2(x i\x2)=ng(\x i-X2\)  inside  the  thick  ply  \x\<aw  and  0  otherwise.  We  will  consider  the  fiber 
volume  fraction  inside  the  ply  c=0.65,  and  two  radial  distribution  functions  g(r)  (step  and  nonstep 
[54]  functions).  The  neglect  of  the  binary  interaction  of  inclusions  for  statistically  homogeneous 
medium  n(Ar)=const.  reduces  the  formula  for  the  effective  elastic  moduli  to  the  analogous  relation 
obtained  by  the  MT  method  which  is  invariant  to  the  g(r).  Assume  the  matrix  is  epoxy  resin 
(A:'^^^=4.27  GPa  and  p*^^^^1.53  GPa)  which  contains  identical  circular  glass  fibers  =50.89  GPa 
and  [1^^=35.QA  GPa).  As  can  be  seen  in  Figure  79,  the  use  of  the  approach  based  on  the 
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Figure  79.  Variation  of  the  Effective  Moduius  L  om  as  a  Function  of  a  Coordinate  Estimated 
by:  MEFM  and  Nonstep  Function  g(r)  (dashed  iine),  MEFM  and  Step  Function  g(r) 
(dot-dashed  iine),  MT  and  Nonstep  grCr)  (dotted  line),  MT  and  Step  Function  gr(r)  (dotted 
line). 

quasicrystalline  approximation  (also  called  the  MT  approach)  leads  to  an  underestimate  of  the 
effective  moduli  L* 2222,  compared  to  the  more  exact  approximation  of  the  MEFM  which  provides 
a  good  comparison  with  experimental  data  for  statistically  homogeneous  media  (see  Figure  76). 
The  distribution  L*(x)=  L* (X2)  was  used  for  estimation  of  statistical  average  stresses  inside  the  ply 
«5>(x)  by  FEA.  The  stress  concentrator  factors  found  at  the  previous  local  evaluations  permit 


one  to  estimate  the  stresses  in  the  inclusions  <c>,(x2)  (see  Figure  80). 


Figure  80.  Variation  of  Statistical  Average  Stresses  Inside  the  Inclusions  «322>i(x2)  (solid  line) 
and  <a72>i(x2)  (dot -dashed  line)  for  the  External  Loading  at  the  Infinity  y=di28j2  and 
aif=8ii5j2,  respectively,  and  a  Step  Function  g. 
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11.  QUANTITATIVE  DESCRIPTION  AND  NUMERICAL  SIMULATION  OE 
RANDOM  MICROSTRUCTURES  OF  COMPOSITES  AND  THEIR 
EFFECTIVE  ELASTIC  MODULI  [59] 

A  digital  image  processing  technique  is  used  for  measurement  of  centroid  coordinates  of 
fibers  with  forthcoming  estimation  of  statistical  parameters  and  functions  describing  the 
stoichiastic  structure  of  a  real  fiber  composite.  Comparative  statistical  analysis  of  the  real  and 
numerically  simulated  structure  was  performed.  Accompanying  known  methods  of  the 
generation  of  random  configurations  by  the  random  shaking  procedure  allows  creation  of  the  most 
homogenized  and  mixed  structures  that  do  not  depend  on  the  initial  protocol  of  particle 
generation.  We  consider  a  linearly  elastic  composite  medium,  which  consists  of  a  homogeneous 
matrix  containing  a  statistically  homogeneous  set  of  ellipsoidal  inclusions.  The  MEFM  (see  for 
reference  [49]),  based  on  the  theory  of  functions  of  random  variables  and  Green's  functions,  is 
used  for  demonstration  of  the  dependence  of  effective  elastic  moduli  of  fiber  composites  on  the 
radial  distribution  functions  estimated  from  the  real  experimental  data  as  well  as  from  the 
ensembles  generated  by  the  method  proposed. 

The  quantitative  description  of  the  microtopology  of  heterogeneous  media  -  such  as 
composite  materials,  porous  and  cracked  solids,  suspensions  and  amorphous  materials  -  is  crucial 
in  the  prediction  of  overall  physical  properties  of  these  materials.  After  many  years  of 
comprehensive  study  by  direct  measurements  and  empirical  relations  that  is  extremely  laborious, 
the  structure  of  microinhomogeneous  materials  is  not  completely  understood.  Computer 
simulation  of  topologically  disordered  structures  by  the  random  packing  of  hard  spherical 
particles  in  2-D  and  3-D  cases  has  a  long  history  originated  in  the  theory  of  liquids.  This  problem 
is  closely  connected  with  the  known  fundamental  problem  of  statistical  physics  -  description  of 
the  behavior  of  the  particle  system  with  the  interaction  potential  of  hard  spheres  (see  for  reference 
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[59]).  The  random  packing  of  spheres  has  been  studied  very  extensively  due  to  the  technological 
importance,  and  the  opportunity  to  model  the  simple  liquid,  concentrated  suspensions,  amorphous 
and  powder  materials. 

Computer  simulation  of  packing  problems  can  be  classified  into  three  groups  of  methods: 
molecular  dynamic,  Monte  Carlo,  and  dense  random  packing.  Much  progress  in  the  theory  of  the 
dense  random  packing  has  been  made  by  the  use  of  two  kinds  of  methods:  the  sequential 
generation  models  and  the  collective  rearrangement  models.  In  the  sequential  model  (see  for 
reference  [59]),  so-called  cluster  growth  model,  a  particle  added  to  the  surface  of  the  particle 
cluster  which  grows  outwards  is  placed  sequentially  to  the  point  closest  to  the  original  such  that 
the  new  particle  established  contact  with  three  existing  spheres  in  the  cluster.  Overlapping  is 
ruled  out  by  checking  the  center-center  distance  of  the  particle.  The  phenomenological  character 
of  the  construction  algorithm  posed  by  the  particle  cluster  from  the  initial  term  containing  three 
particles  leads  to  the  inhomogeneous  and  anisotropic  inclusion  fields  with  a  different  density. 
Moreover,  the  configurations  generated  do  not  demonstrate  the  characteristic  split-second  peak  in 
the  radial  distribution  function  observed  in  the  experimental  packing.  In  the  second  part  of  the 
sequential  generation  models,  called  the  model  of  “rigid  sphere  free-fall  into  a  virtual  box,”  one 
particle  is  dropped  vertically  each  time  from  the  random  point  onto  the  surface  of  an  existing 
particle  cluster  growing  upward.  The  different  densities  were  approached  by  introduction  of  the 
phenomenological  parameter  limiting  the  number  of  rotations  of  a  fallen  sphere  up  to  its  setting  in 
the  structure.  The  effects  of  boundaries  of  the  virtual  box  are  eliminated  by  introducing  the 
conventional  cyclic  boundary  conditions.  The  algorithms  described  belong  to  the  class  of  static 
method  where  the  particles  are  placed  at  a  given  time  step  and  cannot  thereafter  move.  For 
contrast,  the  dynamic  methods  assume  the  reorganization  of  whole  packing  due  to  the  short-  or 
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long-range  interactions  between  particles.  So,  in  collective  rearrangement,  model  N  points 
randomly  distributed  in  a  virtual  box  are  assigned  the  radii  and  the  random  vector  of  moving. 

Each  sphere  is  moved  until  there  are  no  overlaps.  Then  the  radii  are  increased,  and  the  process  is 
repeated  until  any  further  increase  in  radii  or  any  displacement  of  the  spheres  create  overlaps  that 
cannot  be  eliminated  (the  references  on  the  different  versions  of  this  method  can  be  found  in 
[59]).  More  recently,  numerical  simulations  were  performed  to  realize  homogeneous  and 
isotropic  packing  of  spheres  by  various  methods,  for  instance,  assuming  a  hypothetical  sphere 
having  dual  structure  whose  inner  diameter  defines  the  true  density  and  the  outer  one  a  nominal 
density.  An  alternative  approach  eliminating  boundary  effects  of  the  virtual  box  is  based  on  the 
use  of  spherical  boundary  conditions  instead  of  the  periodic  ones.  There,  one  simulates  hard  disks 
(more  exactly,  a  circular  cap  visualized  as  a  contact  lens  on  the  surface  of  an  eyeball)  on  the 
surface  of  the  ordinary  3-D  sphere  and  hard  spheres  on  the  “surface”  of  a  4-D  hypersphere.  The 
advantage  of  this  procedure  is  that  there  is  no  preferential  direction,  and  it  is  impossible  to  pack 
particles  into  perfect  regular  periodic  configurations. 

A  close  random  ensemble  of  spheres  has  been  studied  for  many  years,  and  the  following 
quantitative  parameters  are  well  known.  Packing  densities  for  close  lattice  packing  are  0.9069 
(triangular)  in  the  case  of  disks  packed  into  the  plane  and  0.7405  (fee  or  hep)  in  the  case  of 
spheres  packed  into  1^.  Model  experiments  were  performed  using  steel  balls  of  equal  size 
randomly  packed  into  the  shaking  containers  while  forthcoming  extrapolating  the  measured 
densities  to  eliminate  finite- size  effects  that  provided  the  conventional  value  of  random  close 
packing,  such  as  0.6366.  Random  loose  packing  equal  to  0.60  is  observed  at  the  gentle  rolling  of 
steel  balls  into  the  packing  container  without  shaking.  In  two  dimensions,  the  experimental 
numbers  for  close  and  loose  random  packing  are  estimated  with  less  accuracy,  c=0.8225  and 
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c=0.601  (see  for  reference  [59]).  It  should  be  mentioned  that  in  contradiction  to  the  periodic 
packing,  a  random  close  packing  is  an  ill-defined  problem  which  has  no  unique  theoretical 
definition,  and  its  final  states  are  protocol  dependent  in  both  the  numerical  simulation  and 
experimental  sense  (see  for  details  [50]). 

In  subsection  11.1,  the  quantitative  descriptors  of  the  dispersion  of  fibers  in  unidirectional 
composites  will  be  analyzed  in  order  to  describe  the  pattern  of  fiber  location  as  it  really  exists 
rather  than  as  described  by  some  assumed  model.  Since  the  structure  of  random  packing  is 
strongly  dependent  on  the  procedure  of  their  generation,  we  will  consider  a  few  popular 
algorithms  and  their  combinations  adapted  for  obtaining  most  homogeneous  configurations,  and 
will  compare  the  statistical  parameters  of  configurations  generated  by  the  different  methods.  In 
subsection  11.2  one  will  estimate  the  dependence  of  effective  elastic  properties  of  fiber 
composites  on  the  radial  distribution  functions  estimated  from  the  real  experimental  data  as  well 
as  from  the  ensembles  generated  by  the  method  proposed. 

11.1  STATISTICAL  DESCRIPTION  AND  NUMERICAL  SIMULATION  OE 
RANDOM  STRUCTURE  COMPOSITES 

The  widely  used  informative  function  describing  the  point  distribution  is  a  second-order 
intensity  function  K( r)  defined  as  the  number  of  further  points  expected  to  locate  within  a  distance 
r  of  an  arbitrary  point  dividing  on  the  number  of  points  per  unit  area  n.  Since  points  lying  outside 
the  observation  window  w  are  not  observed,  it  is  known  that  the  estimator  for  K(r)  should  be 
corrected.  The  function  K( r)  is  obtained  by  averaging  over  all  inclusions  at  each  value  of  r.  Due 
to  the  wide  utilization  of  periodic  boundary  conditions  at  the  numerical  simulation  of  the  random 
packing,  an  alternative  toroidal  edge  collection  is  used  in  which  rectangular  regions  w  can  be 
regarded  as  a  torus  so  that  points  on  opposite  edges  are  considered  to  be  closed.  Then  w  is 
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considered  as  part  of  a  grid  of  identical  rectangles,  forming  a  border  of  the  pattern  inside  w,  and 
the  distances  are  measured  from  the  point  in  the  central  rectangle  w  to  the  points  in  the 
surrounding  periodic  rectangles. 

The  radial  distribution  function  (RDF),  g(r),  playing  a  role  similar  to  the  variance  in  a 
classical  analysis  of  random  variables,  is  defined  as  the  radial  distribution  of  the  average  number 
of  sphere  centers  per  unit  area  in  a  spherical  shell  and  can  be  estimated  from  second-order 
intensity  factor  as  g(r)=(2pr)'^dK(r)/dr.  The  RDF  is  related  to  the  derivative  of  the  K(r),  and 
therefore  it  is  more  sensitive  to  changes  in  the  spatial  order  than  is  the  K( r)  function. 

Digital  images  of  transverse  sections  through  the  fiber  composite  material  were  obtained 
by  digitizing  high- resolution  optical  micrographs  using  standard  methods.  Image  processing  was 
then  carried  out  on  the  1024  by  1024  digital  images  using  a  commercially  available  desktop 
software  package  (Adobe  Photoshop  5.5)  in  conjunction  with  a  plug-in  Image  Processing  Tool 
Kit  (IPTK  3.0  by  John  Russ,  1996-1999  Reindeer  Games,  Inc.). 

Since  the  structures  of  random  packing  are  strongly  dependent  on  the  procedure  of  their 
generation,  we  will  consider  a  few  popular  algorithms  and  their  combinations  and  compare  the 
statistical  parameters  of  configurations  generated  by  the  different  methods. 

Poisson  distribution.  If  the  radius  spheres  are  small  enough,  then  their  placement  is 
described  by  the  stationary  (or  homogeneous)  Poisson  point  process.  Although  the  hypothesis  of 
the  Poisson  set  of  centers  of  nonoverlapping  spheres  is  not  fulfilled  for  the  finite  radius  of 
spheres,  it  can  be  sometimes  used  as  a  useful  approximate  description  of  an  observed  pattern.  We 
recall  that  for  Poisson  distribution  (d=2)  K(r)=pP  and  g(r)=l. 

Hard  core  model  (HCM).  In  this  model  (also  called  random  sequential  adsorption),  the 
disks  with  radius  a  are  placed  one  by  one  with  the  center  positions  A=X7,...,a:a?  being  distributed 
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randomly  and  uniformly  over  the  set  of  all  points  in  a  rectangular  region  w.  One  supposes  that 
periodic  boundary  conditions,  that  is,  w  and  X,  are  periodically  replicated  in  all  directions.  If  the 
new  disk  does  not  overlap  already  deposited  disks,  its  position  is  fixed  and  does  not  move 
anymore;  otherwise,  it  is  rejected  and  another  random  center  position  is  generated.  The  process  is 
finished  when  either  a  preassigned  packing  fraction  is  achieved  or  when  no  particles  can  be  added 
anymore  (jamming  limit),  which  occurs  at  a  volume  fraction  c=  0.55  (2-D  case)  or  c=0.38  (3-D 
case).  The  HCM  provides  a  more  realistic  reference  model  than  a  Poisson  point  process,  in  which 
arbitrary  small  distances  between  points  are  allowed.  The  advantage  of  the  HCM  as  a  protocol 
independence  is  sacrificed  in  the  case  of  generation  of  the  binary  or  polydisperse  structures.  The 
geometrical  blocking  effect  and  the  process  irreversibility  leads  to  the  configuration  packing  that 
is  essentially  different  from  corresponding  equilibrium  configuration.  To  prevent  just  this  kind  of 
low- density  jam  from  occurring  in  the  following  model,  we  will  shake  the  disks. 

Hard  core  shaking  model  (HCSM).  This  is  a  sort  of  HCM  generating  increasing  number 
of  inclusions  in  a  virtual  box  w,  accompanied  by  the  shaking  process,  by  giving  each  disk  a  small 
random  displacement  independent  of  its  neighbors'  positions  that  make  it  possible  to  unlock  the 
disks  from  the  jamming  configuration  (taking  place  for  HCM  at  c=0.55),  and  allows  them  to  find 
the  most  homogeneous  and  mixed  arrangement.  Various  algorithms  have  been  devised  to 
simulate  reordering  due  to  shaking  or  vibration  of  dense  packing,  which  reduces  the  volume 
concentration  of  the  high- density  jam  configuration.  The  packing  with  a  wide  concentration  of 
inclusions  has  been  investigated  less  and  necessitates  some  additional  consideration.  In  Figure  8 1 
the  comparative  analysis  of  the  RDF  g(r)  as  a  function  of  r  estimated  by  the  HCM  and  HCSM  is 
presented  for  the  fiber  concentration  c=0.5,  0.55.  The  plot  of  g(r)  of  the  configuration  generated 
by  the  HCM  was  obtained  by  averaging  over  30  realizations;  other  curves  correspond  to  one 
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Figure  81.  The  Radial  Distribution  Functions  g(r)  versus  Relative  Radius  r/a  Estimated  by  the  HCM 
at  c=0.55  (solid  line)  and  c=0.5  (dotted  line),  and  by  HCSM  at  c=0.55  (dot-dashed  line) 
and  c=0.5  (dashed  line). 


realization.  In  so  doing,  the  regular  shaking  procedure  of  the  HCSM  yields  the  fitting  of 
averaging  curves.  As  can  be  seen,  the  HCSM  leads  to  more  long-range  ordering  than  HCM,  and 
at  the  same  c  the  RDF  g( r)  for  a  small  r  is  higher  in  the  case  of  structure  simulation  by  the  HCM 
rather  then  by  HCSM.  It  should  be  mentioned  that  for  a  population  of  finite- size  fibers,  a  value  of 
g{r)  higher  than  1  does  not  necessarily  imply  that  the  fibers  are  clustered.  The  comparison  of  the 
histograms  of  the  average  number  rit  of  inclusions  in  the  rectangular  testing  window  w\={x:  Ix/- 
Xii\<Rt,  1=1,2}  (Rt=3a)  is  presented  in  Figure  82,  from  which  we  notice  that  the  compromise  of 
the  shaking  procedure  with  the  modified  HCM  leads  to  a  more  homogeneous  and  mixed 
arrangement;  the  fraction  p(nt)  of  the  testing  windows  containing  both  the  small  and  large 
numbers  rit  of  inclusions  is  diminished.  Obviously,  the  descriptor  rit  is  a  more  sensitive  measure 
of  the  local  statistical  homogeneity  of  the  configuration  analyzed  than  the  RDF  g(r). 
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Number  of  inclusions 

Figure  82.  Histogram  of  Fractions  p(nt)  of  Testing  Window  Containing  rit  Inclusions  Generated  by 
the  HCM  (solid  curve)  and  by  the  HCSM  (dotted  curve). 

Monte  Carlo  shaking  process.  The  Monte  Carlo  shaking  process  is,  structurally,  the  most 
influential,  and,  computationally,  the  most  intensive  part  of  the  whole  packing  process.  The 
duration  of  this  phase,  which  can  be  measured  in  terms  of  the  number  of  Monte  Carlo  steps  per 
particle,  N^^/N,  can  be  changed  by  variation  of  the  size  of  the  shake-up  window,  the  frequency 
ratio  of  the  alternation  of  the  processes  in  the  generation  of  a  new  particle,  and  the  shaking. 
However,  the  results  of  total  simulation  are  not  related  trivially  to  the  detailed  separate  stage.  In 
Figure  83  we  have  plotted  for  one  plausible  input  parameter  set,  the  volume  fraction  against  the 
length  of  Monte  Carlo  simulation  N^^/N.  The  HCM  allows  one  to  generate  the  random  packing 
more  steeply  than  the  HCSM,  only  at  the  small  particle  concentration  c<0.52.  As  the  packing 
density  grows,  the  percent  of  trial  nonaccepted  generations  dramatically  increases  in  the  vicinity 
of  the  jamming  limit,  and  equals,  for  example,  99.998  percent  at  c=0.5445  and  A=3130  with  CPU 
time  expended  equal  to  five  hours  for  a  PC  with  a  644  MHz  processor.  At  c=0.65  the  advantage 
of  the  HCSM  comprising  in  creating  more  homogeneous  configurations  in  comparison  with  the 
HCM  (see  Figure  82)  degenerates.  Additional  shortcomings  are  the  absence  of  the  testing 
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Figure  83.  Length  of  Monte  Carlo  Simulation  t/^^/Nvs.  an  Area  Fraction  c  Appropriate  to  the  HCM 
(solid  curve)  and  to  the  HCSM  (dotted  curve). 

windows,  with  a  small  number  of  particles  and  a  finite  empty  volume  with  a  diameter  more  than 
2a.  Thus,  the  HCSM  stumbles  on  an  intrinsic  obstacle  in  the  form  of  the  jamming  limits 
^HCHM_Q  ^rj  more  fuzzy  than  the  jamming  limit  for  the  HCM,  c^‘'^=0.55.  Forthcoming 

expansion  of  the  particle  concentration  is  possible  through  the  utilization  of  the  growth  of  the 
particle  radius  considered  [60]. 

Initially  periodic  shaking  model  (IPSM).  The  periodic  lattice  packing  ?  :Xm=miei+m2e2 
of  centers  of  the  disks  with  the  radius  a  was  chosen  as  an  initial  system.  Here  ei  and  ei  are 
linearly  independent  vectors.  The  rearrangement  of  initially  periodic  structure  is  conducted  by 
the  shaking  procedure  described  in  [59].  As  can  be  seen  in  Figure  84,  where  g(r)  is  estimated  by 
IPSM  with  initial  triangular  packing  of  inclusions  with  different  numbers  of  global  shaking 
presented  at  c=0.65,  the  evaluation  of  RDF  becomes  stable  if  the  numbers  of  global  shaking  are 
more  than  50.  In  so  doing,  CPU  time  expended  for  one  global  shaking  equals  10^  sec.  at  N=3700, 
Rt=3.1a.  Comparison  of  g(r),  estimated  by  IPSM  with  initial  triangular  and  square  packing  of 
inclusions  at  100  global  shakings  as  well  as  by  the  HCSM,  was  estimated  also  at  c=0.65,  from 
which  there  is  reason  to  believe  that  a  sufficiently  intensive  shaking  process  eliminates  the 
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Figure  84.  The  Radial  Distribution  Functions  g(r)  versus  Relative  Radius  r/a  Estimated  by  the 
IPSM  at  c=0.65,  N=3700:  t50Shaking  (solid  curve),  tOO  Shaking  (dashed  curve),  30 
Shaking  (dot-dashed  curve),  10  Shaking  (dotted  curve). 

sensitivity  of  results  obtained  on  the  concrete  algorithm  of  the  simulation  of  the  random  inclusion 
ensemble,  as  well  as  on  the  initial  arrangement  of  inclusions.  Estimations  of  the  real g(r)  are 
presented  in  Figure  85. 

Collective  rearrangement  model  ( CRM).  The  CRM  proposed  irv  [60]  leads  to  close 
packing.  Because  our  goal  is  different  and  consists  of  the  simulation  of  well- mixed  random 
structures  in  a  complete  range  of  the  inclusion  concentration,  some  corrections  of  known 
algorithms  summarized  above  were  carried  out.  The  modification  lies  in  the  fact  that  after  a  few 
steps  of  the  estimation  of  the  new  velocities  of  colliding  inclusions,  both  procedures  in  the 
generation  of  new  inclusions  in  low-density  testing  windows  and  of  inclusion  shaking  described 
m  [59]  are  realized.  It  should  be  mentioned  that  modified  CRM  is  not  as  optimal  as  an  original 
CRM  in  the  sense  of  modeling  of  close-packing  configuration,  simply  because  an  added 
procedure  of  random  shaking  is  just  focused  on  the  destruction  of  dense,  locked  local 
configurations  in  some  testing  windows  and  on  the  generation  of  most  homogeneous  structures. 
On  the  other  hand,  the  support  of  the  known  simulation  protocol  by  a  shaking  procedure  has  some 
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Figure  85.  The  Radial  Distribution  Functions  g(r)  versus  Relative  Radius  r/a:  Single  Sample  546- 
0181  (dashed  curve),  Averaged  Over  10  Samples  of  546-0161  (dotted  curve),  at  c=0.65, 
Averaged  Over  10  Samples  of  546-0181  (dot-dashed  curve).  Averaged  Over  Eight 
Materials  Produced  by  the  Different  Technological  Regime  (solid  curve). 

additional  benefits.  So,  in  the  HCM,  CRM,  and  sequential  generation  model,  the  parameters  and 
functions  are  calculated  only  from  one  simulation  of  the  generated  configuration.  Such  data 
should  be  considered  as  just  a  single  realization  of  a  random  generation  process.  To  provide 
statistically  more  reliable  results,  it  would  be  necessary  to  average  several  realizations.  However, 
this  repeating  procedure  is  not  necessary  in  the  protocols  accompanied  by  a  shaking  procedure 
because  a  configuration  generated  by  a  few  global  shakings  can  be  regarded  as  a  separate 
realization.  Experiments  and  simulations  give  a  suggestion  of  a  transition  between  random  and 
ordered  configurations  in  the  vicinity  of  a  density  c=0.8,  where  the  density  increases  much  more 
slowly  beyond  this  point.  The  understanding  of  this  transition  is  more  obvious  after  quantitative 
analysis,  shown  in  Figures  86  and  87,  of  the  RDFs  estimated  by  the  modified  CRM  for  different 
disk  concentrations  of  c=0.60-0.75  and  c=0.75-0.90,  respectively.  As  can  be  seen  for  large  disk 
concentration,  the  plot  of  the  RDF  has  the  characteristic  split-second  peak  observed  in  the 
experimental  packing,  demonstrating  the  availability  of  large  clusters  with  the  triangular  disk 
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Figure  86.  Radial  Distribution  Functions  g(r)  versus  Relative  Radius  r/a  Estimated  by  the  Modified 
CRM  at  c=0.6(7 (dashed  curve),  c=0.65 (dot-dashed  curve),  c=0.70  (dotted  curve),  c=0.75 
(solid  curve);  Af  grows  from  799  to  811. 


Figure  87.  Radial  Distribution  Functions  g(r)  versus  Relative  Radius  r/a  Estimated  by  the  Modified 
CRM  at  N=811  and  c=0.d0  (solid  curve),  c=0.S5  (dotted  curve),  c=0.80  (dot-dashed 
curve),  c=0. 75  (dashed  curve). 


packing.  For  analysis  of  peculiarity  of  the  plots  g( r)~  r  in  Figure  86,  we  will  compare  the 
coordinates  of  their  peaks  with  the  analogous  peaks  for  close  triangular  and  square  packing.  A 
distinguishing  characteristic  of  the  last  regular  packing  is  the  availability  of  the  fixed  peaks 
corresponding  to  r=2a,4a,6a,...  and  floating  subpeaks,  the  location  of  which  depends  on  the 
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specific  stracture  of  the  unit  cell.  So,  for  the  triangular  and  square  packing,  the  coordinates  of  the 
subpeaks  are  r=3.46a,5.32a,...  and  r=2. 83a, 4-44a, 6.3a...,  respectively.  Imperfections,  both  in 
the  identification  of  the  ball  centers  and  in  the  lattice  itself,  broaden  the  peaks  and  raise  the 
heights  of  the  intervening  values.  Thus,  one  can  presuppose  that  the  second  subpeak  in  Figure  85 
at  r=3.47a  is  caused  by  the  influence  of  local  ordering  in  the  form  of  clusters  with  the  triangular 
structure. 

In  Figure  88  we  compared  the  RDF  estimated  from  experimental  data,  numerical 
simulation  by  the  CRM,  as  well  as  the  RDF  proposed  in  [54]  and  having  analytical  representation 
which  takes  into  account  a  neighboring  order  in  the  distribution  of  the  inclusions.  The  so-called 
well- stirred  approximation  of  the  RDF  differs  from  the  RDF  for  Poisson  distribution  by  the 
availability  of  included  volume,  with  the  center  Xi  where  g( r)=0.  Figure  88  shows  a  good  fit  of 
the  RDFs  estimated  from  experimental  data  and  from  numerical  simulation  by  the  modified  CRM 
and  their  large  dissimilarity  from  analytical  representation  of  g( r). 


Figure  88.  Radial  Distribution  Functions  g(r)  versus  Relative  Radius  r/a  Estimated  by  the 

Numerical  Simulation  (solid  curve),  from  Experimental  Data  (dotted  curve),  by  the 
Analytical  Approximation  (547  (dot -dashed  curve). 
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11.2  EFFECTIVE  ELASTIC  PROPERTIES 

For  a  detailed  discussion  and  numerous  references  of  MEFM  and  related  methods,  readers 
are  referred  io  [49].  Just  for  demonstration  of  the  comparison  of  available  experimental  data  with 
the  predicting  opportunity  of  the  method  proposed,  we  will  consider  the  estimation  of  the 
effective  elastic  moduli  L*.  Assume  that  the  matrix  is  epoxy  resin  {k^^^=A.21  GPa  and 
GPa)  that  contains  identical  circular  glass  fibers  (^^^^=50.89  GPa  and  /r^^^=35.04  GPa).  As  can  be 
seen  from  Figure  89,  the  use  of  the  approach  based  on  the  quasicrystalline  approximation  (MT 
approach)  leads  to  an  underestimation  of  the  effective  shear  modulus  by  1.85  times  for  c=0.7, 
compared  to  the  experimental  data  as  well  as  to  the  more  exact  approximation  of  the  MEFM, 
which  provides  a  good  comparison  with  experimental  data  [55],  In  so  doing,  the  best  fitting  is 
provided  by  the  RDF  simulated  by  the  modified  CRM. 


0.0  0.2  0.4  0.6 

Fiber  volume  fraction 

Figure  89.  Variation  of  the  Effective  Shear  Modulus  p*  as  a  Function  of  a  Concentration  of  the 
Inclusions  c.  Experimental  data  o  and  curves  calculated  by  the  MEFM  and  g(r)  [54] 
(solid  line),  by  MEFM  with  the  RDF  simulated  by  the  modified  CRM  (dot -dashed  line), 
by  the  MEFM  with  well  steered  g(r)  (dashed  curve),  and  by  the  Mori-Tanaka  method 
(dotted  line). 

Let  us  now  demonstrate  the  application  of  the  theoretical  results  by  considering  an 
isotropic  composite  made  of  an  incompressible  isotropic  matrix,  filled  with  rigid  circle 
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inclusions  of  one  size.  This  example  was  chosen  deliberately  because  it  provides  the  maximum 
difference  of  predictions  of  effective  elastic  response  estimated  by  various  methods.  In  Figure 
90  the  most  advanced  micromechanical  model  is  analyzed  for  the  different  RDFs.  As  can  be 
seen,  the  effective  shear  moduli  can  differ  by  a  factor  of  two  or  more  for  the  different  RDFs.  In 
so  doing,  the  RDF  simulated  by  the  modified  CRM  provides  the  estimations  of  that  are 

very  close  to  those  obtained  by  the  real  RDF  at  c=0.65.  It  is  interesting  to  note  that  all  RDFs 
lead  to  the  infinite  values  of  p*  for  large  values  of  c,  but  the  simulated  RDF  provides  a  limiting 
upper  value  of  c=0.72. 


Figure  90.  Variation  of  the  Relative  Effective  Shear  Modulus  pVp*”’  as  a  Function  of  a 

Concentration  of  the  Inclusions  c  Estimated  by  the  MEFM  with  g(r)  [54]  (dot -dashed 
line),  by  the  MEFM  with  the  RDF  Simulated  by  the  Modified  CRM  (solid  curve),  by  the 
MEFM  with  Well  Steered  g(r)  (dashed  curve),  and  by  the  MEFM  with  Experimentally 
Estimated  RDF  o. 


11.3  SUMMARY  AND  RECOMMENDATIONS 

It  should  be  mentioned  that  for  computer- simulated  configurations,  the  particle 
reorganization  induced  by  shaking  is  subjected  only  to  geometrical  constraints,  whereas  for  real 
structures,  the  packing  is  far  more  complicated  and  caused  by  the  elastic,  hydrodynamic,  and 
cohesive  forces  and  others.  Our  simulation  technique  is  able  to  isolate  the  fundamental 
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geometrical  constraints  from  other  physics- mechanical  and  chemical  effects,  and  the  results 
provide  a  valuable  benchmark  for  evaluating  sophisticated  packing  schemes  used  to  model  real 
composite  materials.  It  was  shown  that  a  sufficiently  intensive  shaking  process  leads  to  the 
stabilization  of  statistical  distribution  of  the  structure  simulated  that  is  most  homogeneous,  mixed 
and  protocol  independent  (in  the  sense  that  the  statistical  parameters  estimated  do  not  depend  on 
the  basic  simulated  algorithm  such  as  HCM,  CRM,  or  other).  The  arguments  justified  in  the  last 
statement  are  plausible  rather  than  rigorous  and  require  additional  investigation. 

It  is  known  that  taking  only  one  point  probability  density  (volume  fraction)  into  account 
can  provide  just  a  rough  estimation  of  bounds  of  effective  properties  and  statistical  averages  of 
stresses  in  the  constitutive  equations  of  composite  materials.  More  informative  characteristics  of 
the  point  set  are  obtained  using  statistical  second-order  quantities  (such  as  two -point  probability 
density,  second -order  intensity  function,  and  nearest  neighbor  distribution)  which  examine  the 
association  of  a  point  relative  to  other  points.  Few  contributions  have  paid  attention  to  the 
application  of  these  statistical  distributions  for  generation  of  concrete  realization  of  locations  of  a 
final  number  of  interacting  inclusions  with  their  forthcoming  elastic  analysis.  More  rigorous 
estimation  of  statistical  average  of  stress  fields  in  the  constituents  and,  therefore,  of  effective 
elastic  moduli  are  based  on  the  statistical  averaging  of  random  integral  equations  involved  with 
the  infinite  number  of  inclusions  whose  configurations  are  described  by  the  statistical  second - 
order  functions  (see  for  reference  [49,50]).  In  particular,  in  this  study  we  demonstrated  the 
strong  dependence  of  effective  moduli  on  the  concrete  form  of  radial  distribution  function  and 
demonstrated  strong  differences  between  apparently  similar  distributions. 

It  should  be  mentioned,  however,  that  estimation  of  effective  elastic  moduli  is  a  linear 
problem  with  respect  to  the  stress  field  analyzed,  which  is  less  sensitive  to  the  local  stress 
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distribution  than  nonlinear  micromechanical  problems  of  elastoplastic  deformation,  fracture,  and 
fatigue  of  composite  materials.  Therefore,  the  author  [49]  estimated  the  second  moment  of 
stresses  averaged  over  the  volume  of  the  constitutives  by  the  use  of  the  RDF  with  their 
application  for  the  analysis  of  a  wide  class  of  nonlinear  problems.  These  estimations  of  second 
moments  of  stress  are  defined  by  both  the  random  stress  fluctuations  in  the  components  and  the 
inhomogeneity  of  the  stress  fields  in  the  constituents  that  cannot  be  separated  in  the  framework 
of  the  method  proposed.  There  are  a  few  models  (see  for  reference  [59])  based  on  the  idea  that 
at  high  inclusion  concentration,  the  effective  properties  are  dominated  by  the  interaction  forces 
between  neighboring  particles  that  are  proportional  to  the  ratio  of  the  mean  gap  between 
neighboring  particles  to  particle  diameter.  Obviously,  the  use  of  the  average  value  as  a  mean  gap 
instead  of  random  distribution  leads  to  the  loss  of  statistical  information  about  microtopology  of 
the  composite  and  is  conceptually  questionable,  because  the  estimation  problem  of  the  average  of 
the  output  parameter  (such  as  effective  properties)  by  the  use  of  the  average  of  the  random  input 
parameter  (such  as  nearest  neighbor  distance)  is  essentially  nonlinear.  The  approach  hy  [61] 
based  on  the  elastic  solution  for  two  interacting  inclusions  in  an  infinite  matrix  with  forthcoming 
averaging  is  more  in  perspective,  because  it  is  more  sensitive  to  the  local  configuration  of 
inclusions.  In  so  doing,  the  crude  assumption  hy  [61]  that  two  inclusions  are  subjected  to  the 
external  loading  can  be  relaxed  by  the  use  of  more  accurate  estimations  by  [49]  of  the  correlation 
function  of  the  effective  field  acting  on  two  inclusio  ns  and  fixed  in  the  composite  material. 
However,  more  detailed  realization  of  this  scheme  is  beyond  the  scope  of  the  current  study  and 
will  be  considered  in  forthcoming  publications  by  the  authors. 
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LIST  OF  ACRONYMS 


ADM 

axisymmetric  damage  model 

AE 

acoustic  emission 

B-SAM 

B- spline  analysis  method 

CRM 

collective  rearrangement  model 

ENF 

end  notch  flexure 

E-SEM 

environmental  scanning  electron  microscope 

FE 

finite  element 

FEA 

finite  element  analysis 

FEM 

finite  element  method 

FTM 

Fourier  transform  method 

HCM 

hard  core  model 

HCSM 

hard  core  shaking  model 

EM 

iteration  method 

IPSM 

initially  periodic  shaking  model 

LBC 

ligament  boundary  conditions 

MEFM 

multiparticle  effective  field  method 

MT 

Mori-Tanaka 

NDE 

nondestructive  evaluation 

OBC 

overall  loading/boundary  conditions 

RDE 

radial  distribution  function 

SEM 

scanning  electron  microscope 
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